{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 165,
   "id": "10f89794",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from numpy import * ;\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.linalg\n",
    "import os\n",
    "def load_csv(path):\n",
    "    data_read = pd.read_csv(path)\n",
    "    list = data_read.values.tolist()\n",
    "    data = np.array(list)\n",
    "    print(data.shape)\n",
    "    # print(data)\n",
    "    return data\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "id": "434fe25a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(49, 500)\n",
      "(49, 1)\n",
      "(49, 501)\n",
      "(49, 1)\n",
      "(501, 1)\n"
     ]
    }
   ],
   "source": [
    "AFile = \"randn_data_regression_A.csv\"\n",
    "BFile = \"randn_data_regression_b.csv\"\n",
    "A = load_csv(AFile)\n",
    "B = load_csv(BFile)\n",
    "rows = A.shape[0]\n",
    "cols = A.shape[1]\n",
    "# 拓展矩阵第1列设置为1\n",
    "A=mat(np.c_[np.ones((rows,1)),A])\n",
    "rows = A.shape[0]\n",
    "cols = A.shape[1]\n",
    "B = mat(B)\n",
    "print(A.shape)\n",
    "print(B.shape)\n",
    "#创建theta矩阵\n",
    "theta = np.zeros((cols,1))\n",
    "print(theta.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58092e4c",
   "metadata": {},
   "source": [
    "# 1)\t方法一：解线性方程组计算得到线性回归模型的参数\n",
    "\n",
    "#求解线性方程组求解方法\n",
    "#一 python自带函数求解 \n",
    "#报错 因为必须是 方阵, 也就是行数等于列数\n",
    "#theta = scipy.linalg.solve(A,B)\n",
    "#二 使用自己写的迭代法求解线性方程组的方法求解, 我是之前用matlab写的, 所以需要重新写 也必须是方阵?\n",
    "\n",
    "#AX = b  https://blog.csdn.net/skj1995/article/details/90114611 可逆的判断"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "id": "31357ba3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "矩阵不可逆,请采用梯度下降法\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "before = A.T * A\n",
    "if (np.linalg.det(before)!=0):\n",
    "    theta = before.I*A.T*B \n",
    "    #θ=(x^-1X)^-1X^Ty\n",
    "    print(theta)\n",
    "else:\n",
    "    #可以看到A数据的组数比特征少\n",
    "    print(\"矩阵不可逆,请采用梯度下降法\")\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42e8daba",
   "metadata": {},
   "source": [
    "# 2)\t方法二：利用梯度下降法计算线性回归模型的参数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "ec78dc25",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradientDescent(X, y, theta, learning_rate, num_iters, tol,print_step):\n",
    "    #参数说明: X数据, y预测目标值, theta参数,learinig_rate 学习率  num_iters迭代次数 tol精度(0.01)\n",
    "    m = len(y)\n",
    "    J_history = np.zeros((num_iters, 1))\n",
    "\n",
    "    for i in range(num_iters):\n",
    "        n = len(theta)\n",
    "        theta_temp = theta\n",
    "        for j in range(n):\n",
    "            theta_temp[j] = theta[j] - learning_rate/m*(X*theta-y).T*X[:,j];\n",
    "        theta = theta_temp \n",
    "        cost = computeCost(X, y, theta)\n",
    "        J_history[i] = cost\n",
    "        if(i%print_step==0):\n",
    "            print('第%d次迭代, cost = %f' %(i,cost))\n",
    "        if(cost < tol):\n",
    "            print('迭代训练结束,迭代次数:%d, 偏差值cost=%f'%(i,cost))\n",
    "            return (theta,J_history)\n",
    "    \n",
    "    return (theta,J_history)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "927fdf3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#计算偏差值\n",
    "def computeCost(X,y,theta):\n",
    "    # 这里均方误(欧式距离)差求偏差\n",
    "    m = len(y) #number of training examples\n",
    "    \n",
    "    inner = np.power(((X * theta) - y), 2) # (pred - y)^2\n",
    "    \n",
    "    return np.sum(inner) / (2 * m)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "37a279a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第0次迭代, cost = 243.768989\n",
      "第10次迭代, cost = 44.050490\n",
      "第20次迭代, cost = 10.552179\n",
      "第30次迭代, cost = 2.902350\n",
      "第40次迭代, cost = 0.855970\n",
      "第50次迭代, cost = 0.262709\n",
      "第60次迭代, cost = 0.082717\n",
      "第70次迭代, cost = 0.026514\n",
      "第80次迭代, cost = 0.008612\n",
      "第90次迭代, cost = 0.002826\n",
      "第100次迭代, cost = 0.000935\n",
      "迭代训练结束,迭代次数:100, 偏差值cost=0.000935\n"
     ]
    }
   ],
   "source": [
    "# 这里是求解 方法二\n",
    "\n",
    "# 第一部分: 设定传入参数\n",
    "n = rows\n",
    "m = cols\n",
    "X0 = theta\n",
    "\n",
    "A = A\n",
    "B = B\n",
    "theta = np.zeros((m,1))\n",
    "learning_rate = 0.01\n",
    "step = 2000\n",
    "tol = 0.001\n",
    "print_step = 10\n",
    "\n",
    "#不可逆还可以使用加正则化项来做, 先写梯度下降方法吧\n",
    "\n",
    "# 第二部分: 梯度下降方法\n",
    "(theta,J_hist) = gradientDescent(A,B,theta,learning_rate,step,tol,print_step)\n",
    "#print(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "085f64a5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.04634351e-01]\n",
      " [ 1.62400894e-02]\n",
      " [ 9.29961461e-02]\n",
      " [ 6.73010774e-02]\n",
      " [ 7.30270308e-01]\n",
      " [ 8.61464231e-01]\n",
      " [-3.06918121e-01]\n",
      " [-8.99312104e-03]\n",
      " [ 7.95697883e-01]\n",
      " [ 4.17147141e-01]\n",
      " [-6.30807306e-01]\n",
      " [-9.53032181e-02]\n",
      " [ 2.89267106e-01]\n",
      " [ 6.69853466e-01]\n",
      " [-2.12305784e-02]\n",
      " [ 4.93088308e-01]\n",
      " [ 1.35576311e-01]\n",
      " [-2.74558161e-01]\n",
      " [ 4.94648693e-01]\n",
      " [-4.59965800e-02]\n",
      " [ 6.08435859e-01]\n",
      " [ 4.86572210e-01]\n",
      " [ 8.08605003e-01]\n",
      " [ 7.12989582e-01]\n",
      " [-3.02549774e-01]\n",
      " [-8.31883374e-02]\n",
      " [ 1.55373770e-01]\n",
      " [ 7.16539711e-01]\n",
      " [ 3.44417789e-01]\n",
      " [ 7.08536863e-01]\n",
      " [ 4.95973160e-01]\n",
      " [ 3.82168021e-01]\n",
      " [-3.13244724e-01]\n",
      " [ 1.22889921e+00]\n",
      " [ 1.86272622e-01]\n",
      " [ 2.08047527e-02]\n",
      " [ 5.42738922e-01]\n",
      " [ 4.23724920e-01]\n",
      " [ 8.12797121e-01]\n",
      " [ 6.98345100e-01]\n",
      " [ 9.15665448e-02]\n",
      " [ 2.68957265e-02]\n",
      " [ 2.75499784e-01]\n",
      " [ 9.40905345e-01]\n",
      " [ 4.76747893e-02]\n",
      " [ 7.14117737e-01]\n",
      " [ 6.47808839e-01]\n",
      " [ 8.01784937e-02]\n",
      " [ 8.12837739e-01]\n",
      " [ 7.32994728e-01]\n",
      " [ 1.64845760e-01]\n",
      " [ 2.69518667e-01]\n",
      " [-8.53608197e-02]\n",
      " [-1.07182788e-01]\n",
      " [-3.49213213e-01]\n",
      " [ 4.87644116e-01]\n",
      " [ 1.34918280e-01]\n",
      " [-5.68150296e-01]\n",
      " [-1.08044802e-01]\n",
      " [ 3.13240742e-01]\n",
      " [-6.06036749e-03]\n",
      " [ 3.54098877e-02]\n",
      " [ 1.40778115e-02]\n",
      " [-1.03115533e+00]\n",
      " [-3.71938562e-02]\n",
      " [ 5.07033001e-01]\n",
      " [-3.61691296e-01]\n",
      " [ 4.52137418e-01]\n",
      " [-1.09474509e-01]\n",
      " [-3.08541842e-01]\n",
      " [ 2.19853977e-01]\n",
      " [-7.58446153e-02]\n",
      " [-3.71017815e-01]\n",
      " [ 6.43021844e-01]\n",
      " [-1.15989506e-01]\n",
      " [ 1.54799922e-01]\n",
      " [-2.63589798e-01]\n",
      " [-1.78214010e-02]\n",
      " [ 3.05706059e-03]\n",
      " [-1.08837950e+00]\n",
      " [ 1.41932627e-01]\n",
      " [-1.98569746e-01]\n",
      " [ 2.46818235e-01]\n",
      " [-9.83505269e-01]\n",
      " [ 1.81397238e-01]\n",
      " [-4.52125041e-02]\n",
      " [-8.79209660e-01]\n",
      " [ 8.45912591e-02]\n",
      " [ 3.94977827e-02]\n",
      " [-4.97150558e-01]\n",
      " [-3.29392418e-02]\n",
      " [ 6.55792752e-01]\n",
      " [ 1.16666127e-01]\n",
      " [-4.60025778e-01]\n",
      " [ 7.30466876e-01]\n",
      " [-1.38726548e-01]\n",
      " [ 1.81709528e-01]\n",
      " [ 5.98846017e-01]\n",
      " [-4.70808645e-01]\n",
      " [ 3.46655382e-01]\n",
      " [ 4.87464468e-01]\n",
      " [ 3.84166269e-01]\n",
      " [ 3.38805686e-01]\n",
      " [-4.63190658e-01]\n",
      " [-2.84550667e-01]\n",
      " [-7.03311625e-02]\n",
      " [-1.52102529e-01]\n",
      " [ 3.39532835e-01]\n",
      " [ 1.02754025e-01]\n",
      " [-7.78613451e-01]\n",
      " [-8.37842779e-02]\n",
      " [ 1.41339363e-01]\n",
      " [-4.58866036e-01]\n",
      " [ 7.71401576e-01]\n",
      " [-8.99127039e-02]\n",
      " [-1.05179673e-02]\n",
      " [-3.46055497e-01]\n",
      " [ 4.86428952e-01]\n",
      " [ 5.01717590e-01]\n",
      " [ 3.47079674e-01]\n",
      " [-1.27406111e-01]\n",
      " [-6.46208295e-01]\n",
      " [ 6.30526301e-02]\n",
      " [-1.54216062e-01]\n",
      " [-1.31120460e-01]\n",
      " [ 1.61862682e-01]\n",
      " [ 1.41963195e-01]\n",
      " [ 2.58364864e-02]\n",
      " [ 2.98178330e-01]\n",
      " [-1.14368888e-01]\n",
      " [-7.39423582e-02]\n",
      " [ 3.61372936e-01]\n",
      " [-6.21148031e-01]\n",
      " [-1.42896874e-01]\n",
      " [ 2.00293341e-01]\n",
      " [-9.50661295e-02]\n",
      " [-5.95268862e-01]\n",
      " [ 3.89628399e-01]\n",
      " [-2.49936477e-01]\n",
      " [-1.38372528e-01]\n",
      " [-2.00798585e-01]\n",
      " [ 3.93717369e-01]\n",
      " [-5.66187747e-01]\n",
      " [ 1.51228977e-01]\n",
      " [-2.26093356e-01]\n",
      " [ 8.69645984e-02]\n",
      " [ 5.18164306e-02]\n",
      " [-2.85554332e-01]\n",
      " [ 1.71045696e-03]\n",
      " [-4.70169312e-01]\n",
      " [ 2.78897175e-01]\n",
      " [ 1.29947913e-01]\n",
      " [ 2.69159761e-01]\n",
      " [-5.62661048e-01]\n",
      " [ 5.61631046e-02]\n",
      " [-3.27496153e-01]\n",
      " [-2.23248346e-01]\n",
      " [-1.50360076e-01]\n",
      " [ 2.00396953e-01]\n",
      " [-2.94654803e-01]\n",
      " [-2.88258630e-01]\n",
      " [ 6.56709634e-01]\n",
      " [ 2.34980613e-01]\n",
      " [-2.15025217e-01]\n",
      " [ 5.13363255e-01]\n",
      " [ 2.68059498e-01]\n",
      " [ 2.33117122e-01]\n",
      " [-3.86993473e-01]\n",
      " [ 3.76623045e-01]\n",
      " [ 1.80382346e-01]\n",
      " [ 1.78453935e-01]\n",
      " [-4.46977613e-01]\n",
      " [ 9.78862047e-03]\n",
      " [ 2.48191354e-01]\n",
      " [-8.49003054e-01]\n",
      " [-1.44621586e-01]\n",
      " [-4.21693622e-01]\n",
      " [-4.70446250e-01]\n",
      " [ 4.10842527e-01]\n",
      " [-1.22026214e-02]\n",
      " [-4.47805805e-01]\n",
      " [-2.23903638e-01]\n",
      " [ 3.93203272e-04]\n",
      " [-5.76842699e-02]\n",
      " [ 1.40902551e-01]\n",
      " [-2.48352369e-02]\n",
      " [-3.51501497e-01]\n",
      " [ 8.57273609e-02]\n",
      " [ 4.81231682e-01]\n",
      " [ 1.25456916e-01]\n",
      " [-4.56575476e-01]\n",
      " [-2.38410902e-01]\n",
      " [ 3.56740059e-01]\n",
      " [-8.59969446e-03]\n",
      " [-2.95734176e-01]\n",
      " [-1.65323345e-01]\n",
      " [ 4.44459591e-01]\n",
      " [ 1.00763587e-01]\n",
      " [ 1.56989328e-02]\n",
      " [ 4.80448284e-02]\n",
      " [-4.94013714e-03]\n",
      " [-2.09641960e-01]\n",
      " [ 1.85481925e-01]\n",
      " [ 1.61814820e-01]\n",
      " [ 3.79911619e-01]\n",
      " [ 3.57010530e-01]\n",
      " [-5.19667884e-01]\n",
      " [-1.87864230e-01]\n",
      " [-4.13214494e-01]\n",
      " [-2.36948057e-01]\n",
      " [-1.69142275e-01]\n",
      " [ 5.39500868e-02]\n",
      " [ 3.48567626e-01]\n",
      " [-1.07949070e-01]\n",
      " [-3.40952543e-01]\n",
      " [-9.82199963e-02]\n",
      " [ 4.27958398e-01]\n",
      " [ 4.83080543e-01]\n",
      " [-2.93267097e-01]\n",
      " [ 8.75282875e-01]\n",
      " [-1.05550422e-01]\n",
      " [-3.60789500e-01]\n",
      " [-1.60879504e-01]\n",
      " [-3.49645778e-01]\n",
      " [-4.35868728e-02]\n",
      " [-4.64158833e-03]\n",
      " [ 7.01151764e-02]\n",
      " [-2.72001179e-01]\n",
      " [-2.81776033e-02]\n",
      " [-5.67123293e-02]\n",
      " [ 1.84931864e-01]\n",
      " [ 2.14177438e-01]\n",
      " [-4.79852076e-01]\n",
      " [ 3.83869030e-01]\n",
      " [ 3.72907671e-02]\n",
      " [-4.10472667e-01]\n",
      " [-3.19848127e-01]\n",
      " [ 5.40709905e-01]\n",
      " [-4.11642274e-01]\n",
      " [-7.07310524e-01]\n",
      " [ 7.11357158e-02]\n",
      " [ 1.36667241e-01]\n",
      " [-5.32366850e-01]\n",
      " [-1.58119689e-01]\n",
      " [-5.92383538e-01]\n",
      " [-5.33428567e-01]\n",
      " [-6.69626197e-02]\n",
      " [ 2.48119068e-01]\n",
      " [-8.91014598e-02]\n",
      " [ 8.61884233e-01]\n",
      " [ 2.65713548e-01]\n",
      " [-7.20582109e-02]\n",
      " [ 5.67255836e-01]\n",
      " [ 4.58915257e-01]\n",
      " [ 4.55887660e-01]\n",
      " [-4.24023606e-02]\n",
      " [ 4.63571502e-01]\n",
      " [-4.86654024e-03]\n",
      " [ 1.52172010e-01]\n",
      " [-4.94660490e-01]\n",
      " [-8.24645047e-02]\n",
      " [ 6.05159995e-01]\n",
      " [-6.79423577e-02]\n",
      " [-1.37227468e-01]\n",
      " [ 6.15755187e-02]\n",
      " [ 1.00688858e-01]\n",
      " [ 1.56124070e-01]\n",
      " [ 6.12990114e-01]\n",
      " [ 1.17317078e-01]\n",
      " [-1.32125789e-01]\n",
      " [-8.63741937e-01]\n",
      " [-2.50432839e-01]\n",
      " [ 1.58252347e-01]\n",
      " [ 3.00026970e-01]\n",
      " [ 1.37724694e-01]\n",
      " [ 3.21043321e-01]\n",
      " [-1.45050655e-01]\n",
      " [ 1.53376844e-01]\n",
      " [-1.95093142e-02]\n",
      " [ 1.47057294e-01]\n",
      " [ 3.89336020e-01]\n",
      " [-3.80596470e-01]\n",
      " [-3.47889282e-01]\n",
      " [ 9.60549316e-02]\n",
      " [ 3.25676888e-01]\n",
      " [-1.17761687e-01]\n",
      " [ 9.22634678e-03]\n",
      " [ 3.15820245e-01]\n",
      " [ 3.26091723e-01]\n",
      " [-2.78782121e-01]\n",
      " [ 2.16705425e-01]\n",
      " [-3.88775543e-02]\n",
      " [-5.18925844e-01]\n",
      " [-2.09970441e-01]\n",
      " [ 1.65209647e-01]\n",
      " [ 4.81343820e-01]\n",
      " [-9.78655147e-02]\n",
      " [ 2.79966800e-01]\n",
      " [-6.91278273e-02]\n",
      " [-6.79599875e-01]\n",
      " [ 2.19712906e-02]\n",
      " [-2.20309125e-01]\n",
      " [ 7.59916397e-02]\n",
      " [ 2.64690862e-01]\n",
      " [-1.62676080e-01]\n",
      " [ 1.07948368e-01]\n",
      " [ 3.04245316e-01]\n",
      " [-1.59874338e-01]\n",
      " [-2.53055131e-01]\n",
      " [ 4.81912151e-01]\n",
      " [ 1.38411680e-02]\n",
      " [-8.73315304e-02]\n",
      " [-1.10099712e-01]\n",
      " [ 1.46187339e-02]\n",
      " [-9.33948482e-01]\n",
      " [ 1.62276081e-01]\n",
      " [ 7.03839263e-01]\n",
      " [ 5.07612220e-02]\n",
      " [ 2.21581733e-04]\n",
      " [-3.14128059e-01]\n",
      " [-4.02787871e-01]\n",
      " [-4.03503012e-01]\n",
      " [ 6.92546529e-01]\n",
      " [-3.33853016e-01]\n",
      " [-3.39405501e-01]\n",
      " [-1.11920729e-01]\n",
      " [ 1.20935526e-01]\n",
      " [-4.07381089e-01]\n",
      " [ 1.58759705e-01]\n",
      " [ 1.11627361e-01]\n",
      " [ 4.84527082e-02]\n",
      " [-5.99430751e-01]\n",
      " [ 6.58261493e-02]\n",
      " [ 2.32370849e-01]\n",
      " [ 6.04710292e-01]\n",
      " [ 3.32015499e-01]\n",
      " [ 7.74689077e-02]\n",
      " [ 2.45779879e-02]\n",
      " [-9.93123577e-02]\n",
      " [ 3.77298162e-01]\n",
      " [-4.22139666e-01]\n",
      " [-1.46088383e-01]\n",
      " [-8.53746199e-02]\n",
      " [-1.20177929e-01]\n",
      " [-2.03353933e-02]\n",
      " [ 5.85211918e-02]\n",
      " [ 9.01918469e-02]\n",
      " [-1.96194534e-01]\n",
      " [ 3.57229544e-02]\n",
      " [-8.81137545e-02]\n",
      " [-1.00207523e-01]\n",
      " [-3.27064084e-01]\n",
      " [ 5.27978984e-01]\n",
      " [-3.06398177e-01]\n",
      " [-1.65855902e-02]\n",
      " [ 2.09297486e-01]\n",
      " [-7.35584618e-01]\n",
      " [-1.79172612e-01]\n",
      " [-2.08438846e-01]\n",
      " [-3.13339027e-01]\n",
      " [-5.41771760e-01]\n",
      " [-1.81078618e-01]\n",
      " [ 3.86651188e-01]\n",
      " [ 2.11027384e-01]\n",
      " [-1.23559967e-01]\n",
      " [ 9.07672812e-02]\n",
      " [ 1.61731745e-01]\n",
      " [ 6.38602458e-02]\n",
      " [ 3.46579206e-01]\n",
      " [-1.87152410e-01]\n",
      " [-5.34704223e-01]\n",
      " [-2.13315897e-01]\n",
      " [-1.15059362e-01]\n",
      " [-4.70144709e-01]\n",
      " [ 3.25792586e-01]\n",
      " [ 4.77112036e-01]\n",
      " [ 1.97197153e-02]\n",
      " [-5.85357456e-01]\n",
      " [-5.01176270e-02]\n",
      " [-9.92994988e-02]\n",
      " [-1.66841267e-01]\n",
      " [ 1.48304608e-01]\n",
      " [-3.37086338e-01]\n",
      " [-2.18122430e-01]\n",
      " [-3.16772750e-01]\n",
      " [ 2.42488204e-01]\n",
      " [ 3.32966650e-01]\n",
      " [-1.20469917e-01]\n",
      " [ 3.43699760e-01]\n",
      " [ 2.60526237e-01]\n",
      " [ 9.08775573e-02]\n",
      " [ 2.80555737e-01]\n",
      " [ 1.25800785e-01]\n",
      " [-5.63380843e-01]\n",
      " [ 1.59821943e-01]\n",
      " [ 2.27929176e-01]\n",
      " [-1.54897834e-01]\n",
      " [-2.95155150e-01]\n",
      " [ 8.29662783e-01]\n",
      " [ 4.85891992e-01]\n",
      " [-5.27356794e-01]\n",
      " [-5.54111701e-01]\n",
      " [-1.07039153e-01]\n",
      " [ 1.66998664e-01]\n",
      " [-1.62212017e-01]\n",
      " [ 4.03045313e-01]\n",
      " [-4.08415527e-01]\n",
      " [ 5.41539301e-01]\n",
      " [ 9.51456716e-02]\n",
      " [-2.62308687e-01]\n",
      " [ 2.44284365e-01]\n",
      " [-5.03805075e-02]\n",
      " [-3.36743987e-01]\n",
      " [ 1.03318025e-01]\n",
      " [ 6.79794468e-01]\n",
      " [-2.14762992e-01]\n",
      " [-8.12574255e-01]\n",
      " [-3.65703371e-01]\n",
      " [-1.80948059e-01]\n",
      " [ 2.79648101e-01]\n",
      " [ 4.09179357e-01]\n",
      " [ 4.34059267e-01]\n",
      " [ 1.36783333e-01]\n",
      " [-1.04107518e-01]\n",
      " [-4.59066254e-02]\n",
      " [-3.82013897e-01]\n",
      " [-4.25405714e-01]\n",
      " [ 8.75831573e-01]\n",
      " [ 5.58167434e-01]\n",
      " [-7.98998428e-01]\n",
      " [-3.19876195e-01]\n",
      " [-1.45036502e-01]\n",
      " [ 9.01718819e-01]\n",
      " [ 5.62771498e-01]\n",
      " [ 2.96348344e-01]\n",
      " [-3.37642992e-01]\n",
      " [ 4.35558914e-01]\n",
      " [-1.94357105e-01]\n",
      " [-1.10936285e-01]\n",
      " [ 2.65957413e-01]\n",
      " [ 7.43377473e-01]\n",
      " [-2.15387791e-01]\n",
      " [ 4.88384519e-01]\n",
      " [ 9.19600695e-03]\n",
      " [-1.17364971e-01]\n",
      " [-1.12805874e-01]\n",
      " [ 2.67341879e-01]\n",
      " [-3.43512600e-01]\n",
      " [ 7.66492010e-02]\n",
      " [ 5.18699992e-02]\n",
      " [ 4.00226563e-02]\n",
      " [ 2.14835296e-02]\n",
      " [-6.11909797e-01]\n",
      " [ 1.77803460e-01]\n",
      " [ 1.11291400e-01]\n",
      " [ 2.17191869e-01]\n",
      " [ 4.58442715e-01]\n",
      " [ 2.72875604e-01]\n",
      " [-7.75204485e-01]\n",
      " [-4.69583987e-01]\n",
      " [-4.51663525e-01]\n",
      " [-2.04753355e-01]\n",
      " [ 1.08976100e-02]\n",
      " [-1.17432544e-01]\n",
      " [ 1.87616499e-01]\n",
      " [-2.99357712e-01]\n",
      " [-4.89245535e-01]\n",
      " [-2.37978945e-01]\n",
      " [-5.69447177e-02]\n",
      " [ 3.77972793e-01]\n",
      " [-2.94777796e-01]\n",
      " [ 8.23967397e-01]\n",
      " [-2.53809683e-01]\n",
      " [-1.29053682e-01]\n",
      " [ 1.60943116e-01]\n",
      " [ 4.51664064e-01]\n",
      " [-3.03808193e-01]\n",
      " [-1.90090202e-01]\n",
      " [ 4.02811254e-01]\n",
      " [ 3.01354835e-01]\n",
      " [ 1.38880453e-01]\n",
      " [-4.16487672e-01]\n",
      " [-6.48195288e-01]\n",
      " [ 2.37887242e-01]\n",
      " [-2.71690847e-01]\n",
      " [-2.48790952e-01]\n",
      " [-2.07387230e-01]\n",
      " [-3.01136433e-01]\n",
      " [-1.97002316e-01]\n",
      " [-2.62106607e-01]\n",
      " [ 2.37732766e-01]\n",
      " [-2.88200618e-01]\n",
      " [-1.59497366e-02]\n",
      " [-2.39807090e-01]\n",
      " [ 2.46935381e-01]\n",
      " [ 2.12575078e-01]\n",
      " [-3.13498123e-01]\n",
      " [ 1.98397440e-01]\n",
      " [ 1.31526891e-01]\n",
      " [-8.33835140e-02]\n",
      " [ 6.14633186e-01]]\n"
     ]
    }
   ],
   "source": [
    "print(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "c85c531a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x21b69993108>]"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaE0lEQVR4nO3dfXRc9X3n8fd3NJLGen60LT8b40DsGIhRCAnZrLdAeCiNk2zaY5om7Da77jZkN+lyTgrJnm2357BLt3lo2Ka0JNDQLps0C0lxKRCIS5eQpAaZAH7CsYKfkW35UbaFZEnz3T/uFR7LEh49jO7MvZ/XOXNm5nfvjL6+9nz08+/+5nfN3RERkXhJRV2AiIhMPYW7iEgMKdxFRGJI4S4iEkMKdxGRGEpHXQBAS0uLL1q0KOoyRERKysaNGw+7e+to24oi3BctWkRHR0fUZYiIlBQz2z3WtgsOy5jZfDN71sy2mtkWM/tc2P6HZrbfzF4ObzfnvOYuM+s0s+1mdsPU/DFERCRf+fTcB4E73P0lM6sFNprZM+G2r7n7l3N3NrNlwBpgOTAH+JGZvcPdh6aycBERGdsFe+7u3uXuL4WPTwLbgLlv85LVwHfdvd/ddwKdwFVTUayIiORnXLNlzGwR8G5gQ9j0WTN71cweNLPGsG0usDfnZfsY5ZeBma01sw4z6+ju7h5/5SIiMqa8w93MaoBHgc+7ew9wH7AEuALoAr4ynh/s7ve7e7u7t7e2jnqyV0REJiivcDezcoJgf9jdvw/g7gfdfcjds8A3OTv0sh+Yn/PyeWGbiIhMk3xmyxjwALDN3b+a096Ws9tHgc3h43XAGjOrNLPFwFLghakrWURELiSf2TLXAJ8ENpnZy2HbF4FbzewKwIFdwO8AuPsWM/sesJVgps3thZops/3ASf7+lTf47Q8spqm6ohA/QkSkJF0w3N39ecBG2fTE27zmbuDuSdSVl52HT/Fnz3Zy84o2hbuISI6SXlumNlMOwMm+gYgrEREpLiUd7nVhuPf0DUZciYhIcSnpcK/NBKNK6rmLiJwrJuGunruISK4SD/dwWOZN9dxFRHKVdLhXpFNkylOc7FfPXUQkV0mHOwS9d425i4icKwbhnqbnTfXcRURylXy412XK6VHPXUTkHCUf7rWZtGbLiIiMUPLhrp67iMj5Sj/cZ6jnLiIyUsmHu2bLiIicr/TDvTJN30CWM4PZqEsRESkaJR/udTO0MqSIyEglH+7D68toZUgRkbNiEO7quYuIjFTy4V6nlSFFRM5T8uGulSFFRM4Xg3BXz11EZKSSD/ezl9pTz11EZFjJh3uNZsuIiJyn5MO9LGXUVKY1W0ZEJEfJhztoZUgRkZFiEe51mXLNlhERyRGLcFfPXUTkXPEJ93713EVEhsUi3OtmlOs6qiIiOWIR7sGwjHruIiLDYhLu5ZzsG8Tdoy5FRKQoxCLc6zLlDGadNweGoi5FRKQoxCLctb6MiMi5LhjuZjbfzJ41s61mtsXMPhe2N5nZM2a2I7xvDNvNzO41s04ze9XMVhb6D3E23DXuLiIC+fXcB4E73H0ZcDVwu5ktA+4E1rv7UmB9+BzgJmBpeFsL3DflVY8wfKm9E/oik4gIkEe4u3uXu78UPj4JbAPmAquBh8LdHgI+Ej5eDfy1B/4ZaDCztqkuPFdjVQUAx04r3EVEYJxj7ma2CHg3sAGY5e5d4aYDwKzw8Vxgb87L9oVtI99rrZl1mFlHd3f3eOs+R2NV0HM/1ntmUu8jIhIXeYe7mdUAjwKfd/ee3G0ezEEc1zxEd7/f3dvdvb21tXU8Lz1PQ9hz17CMiEggr3A3s3KCYH/Y3b8fNh8cHm4J7w+F7fuB+Tkvnxe2FUxdJk1ZytRzFxEJ5TNbxoAHgG3u/tWcTeuA28LHtwGP5bR/Kpw1czVwImf4piDMjIYZ5RzrVc9dRAQgncc+1wCfBDaZ2cth2xeBe4Dvmdmngd3Ab4TbngBuBjqBXuDfTmXBY2moKue4eu4iIkAe4e7uzwM2xuZrR9nfgdsnWde4NVZVaLaMiEgoFt9QheCkqsbcRUQCsQn3xqpyjmvMXUQEiFO4V6vnLiIyLDbh3lBVTv9gljfPaGVIEZHYhPtbSxCo9y4iEqdw1xIEIiLDYhPu9TOCnrtOqoqIxCjcG6uDnrvCXUQkTuGuMXcRkbfEJtwbqoZ77gp3EZHYhHtluoyqijItHiYiQozCHcL1ZdRzFxGJV7g3aAkCEREgZuGunruISCBW4a6eu4hIIFbhrp67iEggZuFezok3BxjKjuta3SIisROrcG+oqsAdet7U0IyIJFuswn14CQINzYhI0sUq3BveWoJAPXcRSbZYhXtzdRDuR0+r5y4iyRarcG+trQTg8Kn+iCsREYlWrMK9uToI9+6TCncRSbZYhXtFOkX9jHL13EUk8WIV7hAMzSjcRSTpYhfuLTUVGpYRkcSLYbhXcviUZsuISLLFM9zVcxeRhItduLfWVnKyf5C+gaGoSxERiUz8wr1G0yFFRGIX7i21wbdUNWNGRJLsguFuZg+a2SEz25zT9odmtt/MXg5vN+dsu8vMOs1su5ndUKjCx9JakwHUcxeRZMun5/5t4MZR2r/m7leEtycAzGwZsAZYHr7mz82sbKqKzcfZnrtmzIhIcl0w3N39OeBonu+3Gviuu/e7+06gE7hqEvWN2/ASBBqWEZEkm8yY+2fN7NVw2KYxbJsL7M3ZZ1/Ydh4zW2tmHWbW0d3dPYkyzlWRTtFQpSUIRCTZJhru9wFLgCuALuAr430Dd7/f3dvdvb21tXWCZYyupaZSY+4ikmgTCnd3P+juQ+6eBb7J2aGX/cD8nF3nhW3TqqWmQj13EUm0CYW7mbXlPP0oMDyTZh2wxswqzWwxsBR4YXIljp+WIBCRpEtfaAcz+w6wCmgxs33AHwCrzOwKwIFdwO8AuPsWM/sesBUYBG5392n/qmhrrYZlRCTZLhju7n7rKM0PvM3+dwN3T6aoyWqpqeRUuARBpnxaZ2KKiBSF2H1DFbQEgYhIPMM9vJZqt06qikhCxTrcD/Uo3EUkmWIZ7rPrg/VlDvb0RVyJiEg0YhnuTVUVVJSleOPEm1GXIiISiViGeyplzK7PcOCEeu4ikkyxDHcIhma6jivcRSSZYhvuc+ozdPVoWEZEkim24T67fgYHTvSRzXrUpYiITLvYhvuchgwDQ86R01pjRkSSJ7bhPrsumA7ZpRkzIpJAsQ33OQ0zAOjSjBkRSaDYhvvwF5m6jqvnLiLJE9twb64OvsjUpW+pikgCxTbczUxz3UUksWIb7gBt+paqiCRU7MNd68uISBLFO9wbZnCwR19kEpHkiXe41wdfZDp8Wuu6i0iyxDzcg7nuGncXkaSJebgHc93f0Fx3EUmYWIf7/MYqAPYeVbiLSLLEOtzrq8qpy6TZffR01KWIiEyrWIc7wMLmavao5y4iCRP7cF/QXMWeI+q5i0iyxD7cFzZVse/YmwwOZaMuRURk2sQ+3Bc0VTGYdS39KyKJEv9wbw5mzOw+0htxJSIi0yf24b6wuRqAPUcV7iKSHLEP99l1GSrKUpoOKSKJEvtwL0sZ85pmsEfDMiKSILEPdwhOqmrMXUSS5ILhbmYPmtkhM9uc09ZkZs+Y2Y7wvjFsNzO718w6zexVM1tZyOLztbCpij1He3HX0r8ikgz59Ny/Ddw4ou1OYL27LwXWh88BbgKWhre1wH1TU+bkLGiu5lT/IMd6B6IuRURkWlww3N39OeDoiObVwEPh44eAj+S0/7UH/hloMLO2Kap1whY2DU+H1ElVEUmGiY65z3L3rvDxAWBW+HgusDdnv31h23nMbK2ZdZhZR3d39wTLyI/muotI0kz6hKoHA9njHsx29/vdvd3d21tbWydbxtta2FxFyuD17lMF/TkiIsViouF+cHi4Jbw/FLbvB+bn7DcvbItUZbqMRc3V/OKgwl1EkmGi4b4OuC18fBvwWE77p8JZM1cDJ3KGbyJ18cwadhw6GXUZIiLTIp+pkN8BfgZcYmb7zOzTwD3A9Wa2A7gufA7wBPA60Al8E/hMQaqegKWzath1pJczg1odUkTiL32hHdz91jE2XTvKvg7cPtmiCmHpzFqGss6uI6d5x6zaqMsRESmoRHxDFYKeO8AOjbuLSAIkJtyXtNZghsbdRSQREhPumfIyFjRVseOQeu4iEn+JCXeApTNr6NSwjIgkQKLC/eKZtbx++JSupyoisZeocF86s4aBIWeXliEQkZhLVriHM2Y6dVJVRGIuUeF+8cxgxsy2LoW7iMRbosK9qiLNktYaNu8/EXUpIiIFlahwB1gxt57NbyjcRSTeEhfu75pbz8Gefg6d7Iu6FBGRgklcuK+YWw+goRkRibXEhfvyOXWYwaZ9PVGXIiJSMIkL9+rKNBe1VLNJPXcRibHEhTsE4+5bdFJVRGIskeG+Ym49XSf6OHyqP+pSREQKIpHh/q7wpKqGZkQkrhIZ7sMnVV/ZezzqUkRECiKR4V6bKeeSWbV07DoWdSkiIgWRyHAHuGpxEy/tOablf0UklhId7r1nhtjyhua7i0j8JDfcFzUB8OKuoxFXIiIy9RIb7jPrMixsruKFnQp3EYmfxIY7wHsWNdGx+xjuHnUpIiJTKtHhftWiJo6ePsMvu3XRbBGJl0SH+3sWB+PuGzQ0IyIxk+hwX9Rcxey6DD/pPBx1KSIiUyrR4W5mrLqklR/vOKz57iISK4kOd4BVl7Rysm+Ql/Ycj7oUEZEpk/hwv+biFtIp49nth6IuRURkyiQ+3Gsz5Vy5sJF/2t4ddSkiIlNmUuFuZrvMbJOZvWxmHWFbk5k9Y2Y7wvvGqSm1cFZdMpNtXT0c7NFFs0UkHqai5/6v3P0Kd28Pn98JrHf3pcD68HlRW3VJKwD/pKEZEYmJQgzLrAYeCh8/BHykAD9jSl06u5Y59Rme3nIw6lJERKbEZMPdgafNbKOZrQ3bZrl7V/j4ADBrtBea2Voz6zCzju7uaMe7zYybV7Tx3I5uTvQORFqLiMhUmGy4f8DdVwI3Abeb2QdzN3qwaMuoC7e4+/3u3u7u7a2trZMsY/JuuXwOA0POD7ceiLoUEZFJm1S4u/v+8P4Q8APgKuCgmbUBhPclMZB9+bx65jfN4PFXuy68s4hIkZtwuJtZtZnVDj8GPgRsBtYBt4W73QY8Ntkip4OZcctlc/hJ52GOnj4TdTkiIpMymZ77LOB5M3sFeAH4B3d/CrgHuN7MdgDXhc9Lwi2XtTGUdZ7arKEZESlt6Ym+0N1fBy4fpf0IcO1kiorKsrY6Lp5ZwyMb9/Kb710QdTkiIhOW+G+o5jIz1rxnPi/tOc62Ll1bVURKl8J9hI9fOY+KdIr/s2FP1KWIiEyYwn2EhqoKbrmsjR/8fD+n+wejLkdEZEIU7qP4xHsXcKp/kHWvvBF1KSIiE6JwH8XKBY28s62OB5/fSTari2eLSOlRuI/CzPjdVUvYcegUT2/VejMiUnoU7mP41RVtLGqu4hvPdhKsoiAiUjoU7mMoSxn/4V8uYdP+E/x4hy6gLSKlReH+Nj62ch5t9Rm++swv1HsXkZKicH8bFekUv3f9O3h573H+XguKiUgJUbhfwL9eOY9lbXX88ZOv0TcwFHU5IiJ5UbhfQFnK+C+3vJP9x9/kged3Rl2OiEheFO55eP+SFm5cPpt71+/g9e5TUZcjInJBCvc8/dHq5VSmU3zhkVcZ0hebRKTIKdzzNLMuwx/82nI6dh/jr36i4RkRKW4K93H42Mq5XPfOmfzPp7bz8t7jUZcjIjImhfs4mBl/8vHLaa2t5DP/e6MuxyciRUvhPk6N1RX8xW9dyeHTZ7j94ZfoH9T0SBEpPgr3CVgxr557PraCn71+hDu+94pWjhSRojPha6gm3cdWzqP7ZD//48nXaKyq4I9WL8fMoi5LRARQuE/K2g9exNHTZ/jL515nYCjL3R9dQVlKAS8i0VO4T4KZcedNl1KZTnHvP3bS0zfAV379CmZUlEVdmogknMJ9ksyM//yhS6jNlPPfn9zGnqM/5f5PtjOnYUbUpYlIgumE6hT59x+8iG99qp1dh3u55X89z9NbDkRdkogkmMJ9Cl37zln83e3X0FafYe3fbOQLj7zC8V7NhReR6adwn2IXz6zhB5+5ht9dtYRHX9rPr3zl//G3L+7RejQiMq0U7gVQkU7x+zdeyuP/8QNc1FLN7z+6iRv+9Dme2NSlkBeRaWHFcPm49vZ27+joiLqMgnB3ntp8gC8/vZ1fdp9mXuMMPnn1Qn6jfT6N1RVRlyciJczMNrp7+6jbFO7TY3AoyzNbD/Ltn+5iw86jVKZTfPjyOfza5XN435Jmysv0nygRGR+Fe5F57UAPD/10N4+9vJ/eM0PUZdJct2wWNy6fzdVLmqnLlEddooiUAIV7keobGOLHOw7z5OYufrT1ID19g6QMls2p472Lm7lqcROXzatndl1GSxuIyHkU7iVgYCjLi7uOsuH1o2zYeYSf7zlO/2AWgMaqcpbNqWNZWx2LW2pY2FzFgqYq5jTM0HIHIgn2duFesG+omtmNwNeBMuBb7n5PoX5WHJSXpXj/khbev6QFgP7BITbvP8GWN3rY+kYPW7t6eOhnuzkTBn7wGmNuwwwWNFfTVpehtbaSlpoKWmvPPm6sqqA2kyatMX2RRClIuJtZGfAN4HpgH/Cima1z962F+HlxVJku48qFTVy5sOmttqGsc6Cnj91HTrPnSC+7j/ay52gve4708lpXD0dOnxlzqmVVRRm1mTS1mXJqM2nqwvvqijSZ8hSZ8jIqy8uCx+kyKsP7TNhWmS6jLGWUlxllKSOdSpEuM9IpI12WIp0K28uCbbn7pswwgqUaUoaGmESmQaF67lcBne7+OoCZfRdYDSjcJ6EsFfTU5zbM4P1Lzt+ezTrHes/Qfaqf7pPB7XjvACf7BjnZF9z3hPfHe8+w52gvp/sH6RsYon8w+9Yw0HQwA4Mg+MPAz31+9hfCub8Uhn9JBPuAYW+91/k/I79fIqPtNmrbKD9l9P3yr2XU1km+p5SWNe+Zz7/7FxdN+fsWKtznAntznu8D3pu7g5mtBdYCLFiwoEBlJEsqZTTXVNJcU8mls8f/+mzWOTOUpW9giL6BLP2Dwf1w+A9mswxlncEhZzDrDGWzDAx50JZ1BoeyYbszMJR9q93dyTq4Q9YdBxhu4+w2D7dls+G9+7ntw88J2/xs20ij/f9l9P3ye/Ho73d+a74/d7LvOXqjlKKWmsqCvG9kq0K6+/3A/RCcUI2qDjkrlTIyqWAoRkRKW6HOsu0H5uc8nxe2iYjINChUuL8ILDWzxWZWAawB1hXoZ4mIyAgFGZZx90Ez+yzwQ4KpkA+6+5ZC/CwRETlfwcbc3f0J4IlCvb+IiIxN32wREYkhhbuISAwp3EVEYkjhLiISQ0WxKqSZdQO7J/jyFuDwFJZTCKpx8oq9Pij+Gou9Pij+GoutvoXu3jrahqII98kws46xlrwsFqpx8oq9Pij+Gou9Pij+Gou9vlwalhERiSGFu4hIDMUh3O+PuoA8qMbJK/b6oPhrLPb6oPhrLPb63lLyY+4iInK+OPTcRURkBIW7iEgMlXS4m9mNZrbdzDrN7M4iqGe+mT1rZlvNbIuZfS5sbzKzZ8xsR3jfWAS1lpnZz83s8fD5YjPbEB7Lvw2Xao6yvgYze8TMXjOzbWb2vmI6jmb2e+Hf8WYz+46ZZaI+hmb2oJkdMrPNOW2jHjML3BvW+qqZrYywxj8J/55fNbMfmFlDzra7whq3m9kNUdSXs+0OM3MzawmfR3IM81Wy4Z5zEe6bgGXArWa2LNqqGATucPdlwNXA7WFNdwLr3X0psD58HrXPAdtynv8x8DV3vxg4Bnw6kqrO+jrwlLtfClxOUGtRHEczmwv8J6Dd3d9FsKz1GqI/ht8GbhzRNtYxuwlYGt7WAvdFWOMzwLvc/TLgF8BdAOFnZw2wPHzNn4ef++muDzObD3wI2JPTHNUxzE9wLcrSuwHvA36Y8/wu4K6o6xpR42PA9cB2oC1sawO2R1zXPIIP+q8AjxNcf/kwkB7t2EZQXz2wk/CEf057URxHzl4juIlg2ezHgRuK4RgCi4DNFzpmwF8Ct46233TXOGLbR4GHw8fnfKYJrg/xvijqAx4h6GTsAlqiPob53Eq2587oF+GeG1Et5zGzRcC7gQ3ALHfvCjcdAGZFVVfoT4EvANnweTNw3N0Hw+dRH8vFQDfwV+HQ0bfMrJoiOY7uvh/4MkEvrgs4AWykuI7hsLGOWbF+fn4beDJ8XBQ1mtlqYL+7vzJiU1HUN5ZSDveiZWY1wKPA5929J3ebB7/iI5t/ama3AIfcfWNUNeQhDawE7nP3dwOnGTEEE+VxDMetVxP8EpoDVDPKf+WLTdT/9i7EzL5EMLT5cNS1DDOzKuCLwH+NupbxKuVwL8qLcJtZOUGwP+zu3w+bD5pZW7i9DTgUVX3ANcCHzWwX8F2CoZmvAw1mNnxlrqiP5T5gn7tvCJ8/QhD2xXIcrwN2unu3uw8A3yc4rsV0DIeNdcyK6vNjZv8GuAX4RPhLCIqjxiUEv8RfCT8z84CXzGx2kdQ3plIO96K7CLeZGfAAsM3dv5qzaR1wW/j4NoKx+Ei4+13uPs/dFxEcs390908AzwIfD3eLusYDwF4zuyRsuhbYSvEcxz3A1WZWFf6dD9dXNMcwx1jHbB3wqXDGx9XAiZzhm2llZjcSDBN+2N17czatA9aYWaWZLSY4cfnCdNbm7pvcfaa7Lwo/M/uAleG/0aI5hqOKetB/kic+biY4u/5L4EtFUM8HCP7b+yrwcni7mWBMez2wA/gR0BR1rWG9q4DHw8cXEXxwOoH/C1RGXNsVQEd4LP8OaCym4wj8N+A1YDPwN0Bl1McQ+A7BOYABghD69FjHjOAk+jfCz84mgpk/UdXYSTB2PfyZ+Yuc/b8U1rgduCmK+kZs38XZE6qRHMN8b1p+QEQkhkp5WEZERMagcBcRiSGFu4hIDCncRURiSOEuIhJDCncRkRhSuIuIxND/BxPmlngxgHCmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(J_hist.size)\n",
    "\n",
    "x_axis = np.linspace(0,150,150)\n",
    "plt.plot(x_axis,J_hist[0:150])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63c0c4a7",
   "metadata": {},
   "source": [
    "# 3)\t在线性回归模型中加入正则化项 ，实现计算过程。 \n",
    "\n",
    "#正则化之后 对于 梯度下降法和 normalequation 正规方程方法算法都有所改变\n",
    "\n",
    "#好处是对于梯度下降 可以避免guonihe\n",
    "\n",
    "#对于normalequation 可以直接求出,不可逆的问题消失了 数学上的证明 比较高深..\n",
    "\n",
    "#下面写求解\n",
    "#先给出正规方程的求解\n",
    "#再修改函数适应正则化\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e371ed1",
   "metadata": {},
   "source": [
    "## 正规方程求解\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "id": "a860e9a3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "251001\n",
      "嘿嘿,加上就可逆了,瞬间求出,看看数量级基本相同参数位置和梯度下降方法求解的是一致的\n",
      "[[ 1.67555831e-01]\n",
      " [ 1.46175195e-02]\n",
      " [ 7.20965336e-02]\n",
      " [ 6.65135545e-02]\n",
      " [ 5.67697138e-01]\n",
      " [ 6.79904898e-01]\n",
      " [-2.34303042e-01]\n",
      " [-2.30868803e-02]\n",
      " [ 6.29242502e-01]\n",
      " [ 3.03330175e-01]\n",
      " [-4.94058931e-01]\n",
      " [-5.41888241e-02]\n",
      " [ 1.92159122e-01]\n",
      " [ 4.99244977e-01]\n",
      " [-2.24109185e-02]\n",
      " [ 3.84054581e-01]\n",
      " [ 1.05943750e-01]\n",
      " [-2.19200806e-01]\n",
      " [ 3.83785603e-01]\n",
      " [-7.07126742e-02]\n",
      " [ 4.70979633e-01]\n",
      " [ 3.41743916e-01]\n",
      " [ 6.38267739e-01]\n",
      " [ 5.29331672e-01]\n",
      " [-2.51965499e-01]\n",
      " [-7.62135864e-02]\n",
      " [ 9.75018126e-02]\n",
      " [ 5.39397898e-01]\n",
      " [ 2.56357689e-01]\n",
      " [ 5.29367622e-01]\n",
      " [ 3.85564877e-01]\n",
      " [ 3.15662723e-01]\n",
      " [-2.02952710e-01]\n",
      " [ 9.59870453e-01]\n",
      " [ 1.58563998e-01]\n",
      " [-1.99937995e-02]\n",
      " [ 4.45363672e-01]\n",
      " [ 3.27571330e-01]\n",
      " [ 6.03582653e-01]\n",
      " [ 5.68254903e-01]\n",
      " [ 5.86329766e-02]\n",
      " [ 4.70873829e-02]\n",
      " [ 2.16539595e-01]\n",
      " [ 7.13162747e-01]\n",
      " [-7.67660422e-04]\n",
      " [ 5.46172125e-01]\n",
      " [ 5.25861277e-01]\n",
      " [ 7.61860113e-02]\n",
      " [ 6.53468050e-01]\n",
      " [ 5.55388398e-01]\n",
      " [ 1.26550907e-01]\n",
      " [ 1.98703951e-01]\n",
      " [-4.50893936e-02]\n",
      " [-7.40597995e-02]\n",
      " [-2.85198725e-01]\n",
      " [ 3.85717861e-01]\n",
      " [ 9.46588427e-02]\n",
      " [-4.62045566e-01]\n",
      " [-1.11580111e-01]\n",
      " [ 2.28363017e-01]\n",
      " [-8.56923730e-03]\n",
      " [ 2.43625711e-02]\n",
      " [ 2.55885731e-02]\n",
      " [-8.13079564e-01]\n",
      " [-4.01266619e-02]\n",
      " [ 4.02722708e-01]\n",
      " [-2.82318404e-01]\n",
      " [ 3.95499622e-01]\n",
      " [-9.19168865e-02]\n",
      " [-2.06025258e-01]\n",
      " [ 1.61240482e-01]\n",
      " [-7.62748825e-02]\n",
      " [-3.01234469e-01]\n",
      " [ 4.98902501e-01]\n",
      " [-8.65037808e-02]\n",
      " [ 1.04089886e-01]\n",
      " [-1.80207951e-01]\n",
      " [-1.94401500e-03]\n",
      " [ 1.19145719e-02]\n",
      " [-8.36076718e-01]\n",
      " [ 1.43313909e-01]\n",
      " [-1.79705778e-01]\n",
      " [ 1.87955416e-01]\n",
      " [-7.74149638e-01]\n",
      " [ 1.28216066e-01]\n",
      " [-5.71776778e-02]\n",
      " [-6.99384721e-01]\n",
      " [ 6.20175394e-02]\n",
      " [ 7.05302483e-02]\n",
      " [-3.42776685e-01]\n",
      " [-3.22252222e-03]\n",
      " [ 5.15815881e-01]\n",
      " [ 8.23537012e-02]\n",
      " [-3.41415351e-01]\n",
      " [ 5.69169714e-01]\n",
      " [-1.03462817e-01]\n",
      " [ 1.65022873e-01]\n",
      " [ 4.76665555e-01]\n",
      " [-4.00817112e-01]\n",
      " [ 2.89289820e-01]\n",
      " [ 3.63965605e-01]\n",
      " [ 3.20282983e-01]\n",
      " [ 2.58613952e-01]\n",
      " [-3.69536263e-01]\n",
      " [-2.00346730e-01]\n",
      " [-5.82566326e-02]\n",
      " [-1.18304808e-01]\n",
      " [ 2.63409436e-01]\n",
      " [ 7.13998617e-02]\n",
      " [-6.18247345e-01]\n",
      " [-5.63505737e-02]\n",
      " [ 9.82152398e-02]\n",
      " [-3.58617652e-01]\n",
      " [ 5.91187553e-01]\n",
      " [-9.17567795e-02]\n",
      " [ 6.82497870e-03]\n",
      " [-2.70384372e-01]\n",
      " [ 3.81086114e-01]\n",
      " [ 4.10007939e-01]\n",
      " [ 2.76261899e-01]\n",
      " [-8.44082875e-02]\n",
      " [-5.34092217e-01]\n",
      " [ 2.07955729e-02]\n",
      " [-1.31879891e-01]\n",
      " [-1.24599893e-01]\n",
      " [ 1.48776215e-01]\n",
      " [ 9.46879224e-02]\n",
      " [ 1.25835052e-02]\n",
      " [ 2.04697936e-01]\n",
      " [-8.38979161e-02]\n",
      " [-3.84704564e-02]\n",
      " [ 2.78656864e-01]\n",
      " [-4.72194620e-01]\n",
      " [-1.25767327e-01]\n",
      " [ 1.47699055e-01]\n",
      " [-9.71567941e-02]\n",
      " [-4.90624927e-01]\n",
      " [ 3.19436868e-01]\n",
      " [-2.03977012e-01]\n",
      " [-9.01506535e-02]\n",
      " [-1.56881963e-01]\n",
      " [ 2.87016368e-01]\n",
      " [-4.37977429e-01]\n",
      " [ 1.66115969e-01]\n",
      " [-1.84247833e-01]\n",
      " [ 6.88292608e-02]\n",
      " [ 2.16766882e-02]\n",
      " [-2.02481970e-01]\n",
      " [-3.88900587e-03]\n",
      " [-3.93385763e-01]\n",
      " [ 2.07680114e-01]\n",
      " [ 9.94774907e-02]\n",
      " [ 2.00092681e-01]\n",
      " [-4.53187533e-01]\n",
      " [ 3.97155170e-02]\n",
      " [-2.23727448e-01]\n",
      " [-1.72492096e-01]\n",
      " [-1.11114692e-01]\n",
      " [ 1.56258619e-01]\n",
      " [-2.51127067e-01]\n",
      " [-2.16315205e-01]\n",
      " [ 5.19237820e-01]\n",
      " [ 1.75828336e-01]\n",
      " [-1.68685748e-01]\n",
      " [ 4.17007179e-01]\n",
      " [ 2.06591550e-01]\n",
      " [ 1.80256573e-01]\n",
      " [-3.05317225e-01]\n",
      " [ 3.23841897e-01]\n",
      " [ 1.44499780e-01]\n",
      " [ 9.66268295e-02]\n",
      " [-3.43733408e-01]\n",
      " [ 2.05439784e-02]\n",
      " [ 2.18210506e-01]\n",
      " [-6.86528532e-01]\n",
      " [-1.04774752e-01]\n",
      " [-3.24789675e-01]\n",
      " [-3.85901013e-01]\n",
      " [ 3.57462803e-01]\n",
      " [-2.24769563e-02]\n",
      " [-3.48028680e-01]\n",
      " [-1.55741462e-01]\n",
      " [ 7.02744828e-03]\n",
      " [-5.63298745e-02]\n",
      " [ 1.24514502e-01]\n",
      " [-4.09226528e-02]\n",
      " [-2.63273124e-01]\n",
      " [ 2.39031808e-02]\n",
      " [ 3.85676353e-01]\n",
      " [ 1.05477143e-01]\n",
      " [-3.51386129e-01]\n",
      " [-1.68080058e-01]\n",
      " [ 3.14188764e-01]\n",
      " [-2.93602574e-02]\n",
      " [-2.27123584e-01]\n",
      " [-1.53860113e-01]\n",
      " [ 3.42160982e-01]\n",
      " [ 1.27203478e-01]\n",
      " [-3.74144431e-04]\n",
      " [ 3.80980442e-02]\n",
      " [ 6.65770933e-03]\n",
      " [-1.53185126e-01]\n",
      " [ 1.39775215e-01]\n",
      " [ 1.30362930e-01]\n",
      " [ 2.88585133e-01]\n",
      " [ 2.58238485e-01]\n",
      " [-3.93245318e-01]\n",
      " [-1.07329917e-01]\n",
      " [-3.30155429e-01]\n",
      " [-2.13134136e-01]\n",
      " [-1.65247487e-01]\n",
      " [ 5.33238875e-02]\n",
      " [ 2.92529837e-01]\n",
      " [-7.67425157e-02]\n",
      " [-3.04235422e-01]\n",
      " [-8.67811110e-02]\n",
      " [ 3.54260170e-01]\n",
      " [ 3.71221554e-01]\n",
      " [-2.76845574e-01]\n",
      " [ 7.21517110e-01]\n",
      " [-1.00572302e-01]\n",
      " [-3.08527110e-01]\n",
      " [-9.50934485e-02]\n",
      " [-2.84588307e-01]\n",
      " [-4.12651291e-02]\n",
      " [ 1.87647835e-02]\n",
      " [ 7.49946990e-02]\n",
      " [-2.39076493e-01]\n",
      " [-2.97497053e-02]\n",
      " [-3.50786521e-02]\n",
      " [ 1.10387998e-01]\n",
      " [ 1.85335422e-01]\n",
      " [-3.85623608e-01]\n",
      " [ 3.11440540e-01]\n",
      " [ 7.39662994e-02]\n",
      " [-3.38566796e-01]\n",
      " [-2.61795649e-01]\n",
      " [ 4.38091076e-01]\n",
      " [-3.13321811e-01]\n",
      " [-5.68887437e-01]\n",
      " [ 7.20879831e-02]\n",
      " [ 1.25821294e-01]\n",
      " [-4.71422174e-01]\n",
      " [-1.19125256e-01]\n",
      " [-4.84157579e-01]\n",
      " [-4.29519884e-01]\n",
      " [-7.07046861e-02]\n",
      " [ 2.05286320e-01]\n",
      " [-8.21598550e-02]\n",
      " [ 6.88578145e-01]\n",
      " [ 2.12900642e-01]\n",
      " [-7.36042402e-02]\n",
      " [ 4.51957193e-01]\n",
      " [ 3.66460510e-01]\n",
      " [ 3.61255658e-01]\n",
      " [-1.24052053e-02]\n",
      " [ 3.50157463e-01]\n",
      " [ 2.29090660e-02]\n",
      " [ 9.33946359e-02]\n",
      " [-3.91221216e-01]\n",
      " [-6.70548330e-02]\n",
      " [ 4.79556725e-01]\n",
      " [-9.21547256e-02]\n",
      " [-1.30651109e-01]\n",
      " [ 5.24797358e-02]\n",
      " [ 9.15773758e-02]\n",
      " [ 1.10238666e-01]\n",
      " [ 5.03803399e-01]\n",
      " [ 8.62663574e-02]\n",
      " [-1.05461291e-01]\n",
      " [-6.76278730e-01]\n",
      " [-1.70249888e-01]\n",
      " [ 1.19113627e-01]\n",
      " [ 2.36147357e-01]\n",
      " [ 9.95541698e-02]\n",
      " [ 2.36136430e-01]\n",
      " [-1.27973915e-01]\n",
      " [ 1.72391755e-01]\n",
      " [ 1.14695277e-02]\n",
      " [ 9.79157639e-02]\n",
      " [ 2.71776514e-01]\n",
      " [-3.35222038e-01]\n",
      " [-2.48481664e-01]\n",
      " [ 1.13902674e-01]\n",
      " [ 2.81897726e-01]\n",
      " [-1.08465031e-01]\n",
      " [ 1.61398151e-02]\n",
      " [ 2.35456901e-01]\n",
      " [ 2.64605266e-01]\n",
      " [-2.07747621e-01]\n",
      " [ 1.86493106e-01]\n",
      " [-8.53104327e-03]\n",
      " [-3.98348245e-01]\n",
      " [-1.67303925e-01]\n",
      " [ 1.18974974e-01]\n",
      " [ 3.78828582e-01]\n",
      " [-1.20065269e-01]\n",
      " [ 2.73572074e-01]\n",
      " [-3.67967665e-02]\n",
      " [-5.71332164e-01]\n",
      " [-1.40023241e-02]\n",
      " [-1.54941513e-01]\n",
      " [ 6.07594719e-02]\n",
      " [ 2.30175777e-01]\n",
      " [-1.57966113e-01]\n",
      " [ 6.65160985e-02]\n",
      " [ 2.85128340e-01]\n",
      " [-1.51641011e-01]\n",
      " [-2.12509735e-01]\n",
      " [ 3.99157561e-01]\n",
      " [ 2.10105193e-02]\n",
      " [-8.70171609e-02]\n",
      " [-6.31815680e-02]\n",
      " [ 1.08399145e-02]\n",
      " [-7.35781732e-01]\n",
      " [ 1.17273064e-01]\n",
      " [ 5.35852671e-01]\n",
      " [ 2.96213641e-02]\n",
      " [-1.84975319e-02]\n",
      " [-2.38753232e-01]\n",
      " [-3.03437279e-01]\n",
      " [-3.39578724e-01]\n",
      " [ 5.77512107e-01]\n",
      " [-2.50131591e-01]\n",
      " [-2.55761458e-01]\n",
      " [-2.85495875e-02]\n",
      " [ 9.72473568e-02]\n",
      " [-3.03620085e-01]\n",
      " [ 1.41200135e-01]\n",
      " [ 1.20239762e-01]\n",
      " [ 6.27652707e-02]\n",
      " [-4.78724380e-01]\n",
      " [ 7.47178713e-02]\n",
      " [ 1.68233500e-01]\n",
      " [ 5.22813012e-01]\n",
      " [ 2.77196638e-01]\n",
      " [ 6.93058504e-02]\n",
      " [ 3.71044932e-02]\n",
      " [-1.09101817e-01]\n",
      " [ 2.64978010e-01]\n",
      " [-3.58777668e-01]\n",
      " [-1.13443850e-01]\n",
      " [-4.04672376e-02]\n",
      " [-1.16245760e-01]\n",
      " [-1.25822956e-03]\n",
      " [ 7.52074238e-02]\n",
      " [ 6.34064785e-02]\n",
      " [-1.49144190e-01]\n",
      " [ 3.63734922e-02]\n",
      " [-7.48006977e-02]\n",
      " [-8.61073522e-02]\n",
      " [-2.41597869e-01]\n",
      " [ 4.52297354e-01]\n",
      " [-2.40278403e-01]\n",
      " [ 2.38398748e-03]\n",
      " [ 1.71210959e-01]\n",
      " [-5.61976985e-01]\n",
      " [-1.22066689e-01]\n",
      " [-1.64363945e-01]\n",
      " [-2.59451976e-01]\n",
      " [-4.10312107e-01]\n",
      " [-9.81426827e-02]\n",
      " [ 2.77361384e-01]\n",
      " [ 1.15983950e-01]\n",
      " [-9.70681333e-02]\n",
      " [ 3.77171932e-02]\n",
      " [ 8.67420977e-02]\n",
      " [ 8.03523211e-02]\n",
      " [ 2.41359805e-01]\n",
      " [-1.59778266e-01]\n",
      " [-4.47081537e-01]\n",
      " [-1.43967190e-01]\n",
      " [-6.40549018e-02]\n",
      " [-3.80933147e-01]\n",
      " [ 2.60646884e-01]\n",
      " [ 3.84245985e-01]\n",
      " [-8.61071994e-03]\n",
      " [-4.52190096e-01]\n",
      " [-6.02650220e-02]\n",
      " [-5.99709627e-02]\n",
      " [-1.67736645e-01]\n",
      " [ 1.28652467e-01]\n",
      " [-2.73772628e-01]\n",
      " [-1.66536755e-01]\n",
      " [-2.52003927e-01]\n",
      " [ 2.13423412e-01]\n",
      " [ 2.89206291e-01]\n",
      " [-8.83764486e-02]\n",
      " [ 2.83018926e-01]\n",
      " [ 2.40824969e-01]\n",
      " [ 5.74640686e-02]\n",
      " [ 2.15545905e-01]\n",
      " [ 1.02737243e-01]\n",
      " [-4.06971043e-01]\n",
      " [ 1.13271823e-01]\n",
      " [ 1.44733771e-01]\n",
      " [-1.36203125e-01]\n",
      " [-2.48258174e-01]\n",
      " [ 6.84929116e-01]\n",
      " [ 3.77996911e-01]\n",
      " [-4.41439535e-01]\n",
      " [-4.12241249e-01]\n",
      " [-9.29910153e-02]\n",
      " [ 1.14104352e-01]\n",
      " [-1.65649466e-01]\n",
      " [ 3.46555660e-01]\n",
      " [-2.98219335e-01]\n",
      " [ 4.20144485e-01]\n",
      " [ 6.18948582e-02]\n",
      " [-2.03981036e-01]\n",
      " [ 2.15574432e-01]\n",
      " [-3.19679687e-02]\n",
      " [-2.50915243e-01]\n",
      " [ 9.24847661e-02]\n",
      " [ 5.45955137e-01]\n",
      " [-1.97488786e-01]\n",
      " [-6.36270577e-01]\n",
      " [-2.73323671e-01]\n",
      " [-1.86356749e-01]\n",
      " [ 2.25975185e-01]\n",
      " [ 3.45193089e-01]\n",
      " [ 3.43184615e-01]\n",
      " [ 1.31381579e-01]\n",
      " [-8.71167742e-02]\n",
      " [-2.76618531e-02]\n",
      " [-3.08793303e-01]\n",
      " [-3.86717889e-01]\n",
      " [ 6.81786282e-01]\n",
      " [ 4.59459534e-01]\n",
      " [-6.48466570e-01]\n",
      " [-2.61614887e-01]\n",
      " [-9.20600194e-02]\n",
      " [ 7.22729465e-01]\n",
      " [ 4.69100353e-01]\n",
      " [ 2.53865622e-01]\n",
      " [-2.45195659e-01]\n",
      " [ 3.27477441e-01]\n",
      " [-1.60341123e-01]\n",
      " [-9.50086265e-02]\n",
      " [ 2.15625310e-01]\n",
      " [ 6.02941811e-01]\n",
      " [-1.31912933e-01]\n",
      " [ 3.93511465e-01]\n",
      " [-1.82568243e-02]\n",
      " [-1.28614380e-01]\n",
      " [-9.85892125e-02]\n",
      " [ 1.94376691e-01]\n",
      " [-2.61246203e-01]\n",
      " [ 1.01232939e-01]\n",
      " [ 4.71645534e-02]\n",
      " [ 6.62601344e-03]\n",
      " [ 2.75725271e-02]\n",
      " [-5.26951760e-01]\n",
      " [ 1.11737531e-01]\n",
      " [ 1.13535240e-01]\n",
      " [ 2.34957721e-01]\n",
      " [ 3.62725219e-01]\n",
      " [ 2.27208909e-01]\n",
      " [-6.57490403e-01]\n",
      " [-4.44949465e-01]\n",
      " [-3.90317154e-01]\n",
      " [-1.23479038e-01]\n",
      " [-8.99600866e-03]\n",
      " [-9.15725848e-02]\n",
      " [ 1.56820478e-01]\n",
      " [-2.94442959e-01]\n",
      " [-3.78863778e-01]\n",
      " [-2.02538467e-01]\n",
      " [-3.99343634e-02]\n",
      " [ 3.21427221e-01]\n",
      " [-2.44195127e-01]\n",
      " [ 6.35825692e-01]\n",
      " [-2.10892486e-01]\n",
      " [-7.44096033e-02]\n",
      " [ 1.55331536e-01]\n",
      " [ 3.22572548e-01]\n",
      " [-2.26323975e-01]\n",
      " [-1.77575348e-01]\n",
      " [ 3.51297388e-01]\n",
      " [ 2.27566278e-01]\n",
      " [ 1.32502746e-01]\n",
      " [-3.04307070e-01]\n",
      " [-5.31918265e-01]\n",
      " [ 2.29621469e-01]\n",
      " [-1.87366660e-01]\n",
      " [-2.04615820e-01]\n",
      " [-1.68671043e-01]\n",
      " [-2.37055758e-01]\n",
      " [-1.68483636e-01]\n",
      " [-2.25065051e-01]\n",
      " [ 1.98779145e-01]\n",
      " [-2.11789679e-01]\n",
      " [ 9.44296040e-03]\n",
      " [-1.78457456e-01]\n",
      " [ 2.27201265e-01]\n",
      " [ 1.77004656e-01]\n",
      " [-2.23322941e-01]\n",
      " [ 1.91546237e-01]\n",
      " [ 1.18904924e-01]\n",
      " [-8.04092873e-02]\n",
      " [ 5.27173567e-01]]\n"
     ]
    }
   ],
   "source": [
    "Regularization_Parameter = 100\n",
    "\n",
    "\n",
    "before = A.T * A+Regularization_Parameter*np.eye(cols)\n",
    "print(before.size)\n",
    "if (np.linalg.det(before)!=0):\n",
    "    print('嘿嘿,加上就可逆了,瞬间求出,看看数量级基本相同参数位置和梯度下降方法求解的是一致的')\n",
    "    theta = before.I*A.T*B \n",
    "    #θ=(x^-1X)^-1X^Ty\n",
    "    print(theta)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "961e0bcf",
   "metadata": {},
   "source": [
    "## 梯度下降求解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "f046216f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#计算偏差值\n",
    "def computeCostWithRegularization(X,y,theta,Regularization_Parameter):\n",
    "    # 加正则化项的计算cost\n",
    "    m = len(y) #number of training examples\n",
    "    \n",
    "    inner = np.power(((X * theta) - y), 2)+ Regularization_Parameter*sum(np.power(theta, 2))\n",
    "    \n",
    "    return np.sum(inner) / (2 * m)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "id": "c5051c74",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradientDescentWithRegularization(X, y, theta, learning_rate, num_iters, tol,print_step,Regularization_Parameter):\n",
    "    #参数说明: X数据, y预测目标值, theta参数,learinig_rate 学习率  num_iters迭代次数 tol精度(0.01)\n",
    "    m = len(y)\n",
    "    J_history = np.zeros((num_iters, 1))\n",
    "\n",
    "    for i in range(num_iters):\n",
    "        n = len(theta)\n",
    "        theta_temp = theta\n",
    "        for j in range(n):\n",
    "            theta_temp[j] = theta[j] - learning_rate/m*(X*theta-y).T*X[:,j];\n",
    "        theta = theta_temp \n",
    "        cost = computeCostWithRegularization(X, y, theta,Regularization_Parameter)\n",
    "        J_history[i] = cost\n",
    "        if(i%print_step==0):\n",
    "            print('第%d次迭代, cost = %f' %(i,cost))\n",
    "        if((i>1) and (abs(J_history[i-1]-cost) < tol)):\n",
    "            print(J_history[i-1]-cost)\n",
    "            print(tol)\n",
    "            print('迭代训练结束,迭代次数:%d, 偏差值cost=%f'%(i,cost))\n",
    "            return (theta,J_history)\n",
    "    \n",
    "    return (theta,J_history)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "id": "6c55f6e8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第0次迭代, cost = 270.030335\n",
      "第10次迭代, cost = 1301.201020\n",
      "第20次迭代, cost = 2273.823487\n",
      "第30次迭代, cost = 2820.044454\n",
      "第40次迭代, cost = 3111.455585\n",
      "第50次迭代, cost = 3266.868856\n",
      "第60次迭代, cost = 3350.442033\n",
      "第70次迭代, cost = 3395.801613\n",
      "第80次迭代, cost = 3420.633711\n",
      "第90次迭代, cost = 3434.332484\n",
      "第100次迭代, cost = 3441.940282\n",
      "第110次迭代, cost = 3446.190101\n",
      "第120次迭代, cost = 3448.576202\n",
      "第130次迭代, cost = 3449.921854\n",
      "第140次迭代, cost = 3450.683693\n",
      "第150次迭代, cost = 3451.116481\n",
      "第160次迭代, cost = 3451.363083\n",
      "第170次迭代, cost = 3451.503972\n",
      "第180次迭代, cost = 3451.584658\n",
      "第190次迭代, cost = 3451.630965\n",
      "第200次迭代, cost = 3451.657592\n",
      "第210次迭代, cost = 3451.672930\n",
      "[-0.00094901]\n",
      "0.001\n",
      "迭代训练结束,迭代次数:214, 偏差值cost=3451.677060\n"
     ]
    }
   ],
   "source": [
    "\n",
    "A = A\n",
    "B = B\n",
    "theta = np.zeros((m,1))\n",
    "learning_rate = 0.01\n",
    "step = 2000\n",
    "tol = 0.001\n",
    "print_step = 10\n",
    "Regularization_Parameter = 100\n",
    "(theta,J_hist) = gradientDescentWithRegularization(A,B,theta,learning_rate,step,tol,print_step,Regularization_Parameter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "fd235696",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.04669171e-01]\n",
      " [ 1.60765702e-02]\n",
      " [ 9.30660330e-02]\n",
      " [ 6.70146276e-02]\n",
      " [ 7.30890709e-01]\n",
      " [ 8.61995685e-01]\n",
      " [-3.07409945e-01]\n",
      " [-8.82881026e-03]\n",
      " [ 7.95877364e-01]\n",
      " [ 4.18309472e-01]\n",
      " [-6.31432282e-01]\n",
      " [-9.57050955e-02]\n",
      " [ 2.90687553e-01]\n",
      " [ 6.71064809e-01]\n",
      " [-2.11533694e-02]\n",
      " [ 4.93908726e-01]\n",
      " [ 1.35924049e-01]\n",
      " [-2.74997037e-01]\n",
      " [ 4.95204467e-01]\n",
      " [-4.48403787e-02]\n",
      " [ 6.08958458e-01]\n",
      " [ 4.88074430e-01]\n",
      " [ 8.09722526e-01]\n",
      " [ 7.14447558e-01]\n",
      " [-3.02394013e-01]\n",
      " [-8.29142063e-02]\n",
      " [ 1.56319169e-01]\n",
      " [ 7.17370375e-01]\n",
      " [ 3.45037689e-01]\n",
      " [ 7.10569725e-01]\n",
      " [ 4.96844394e-01]\n",
      " [ 3.81603923e-01]\n",
      " [-3.15128550e-01]\n",
      " [ 1.23027981e+00]\n",
      " [ 1.86205390e-01]\n",
      " [ 2.23721112e-02]\n",
      " [ 5.42321318e-01]\n",
      " [ 4.24468717e-01]\n",
      " [ 8.14640149e-01]\n",
      " [ 6.98320002e-01]\n",
      " [ 9.18819547e-02]\n",
      " [ 2.64530245e-02]\n",
      " [ 2.75271431e-01]\n",
      " [ 9.42916746e-01]\n",
      " [ 4.88885922e-02]\n",
      " [ 7.15306003e-01]\n",
      " [ 6.47160595e-01]\n",
      " [ 7.94081671e-02]\n",
      " [ 8.13818077e-01]\n",
      " [ 7.33945314e-01]\n",
      " [ 1.65132536e-01]\n",
      " [ 2.69857437e-01]\n",
      " [-8.62908409e-02]\n",
      " [-1.07511182e-01]\n",
      " [-3.50143055e-01]\n",
      " [ 4.87653987e-01]\n",
      " [ 1.35287152e-01]\n",
      " [-5.68183369e-01]\n",
      " [-1.07626792e-01]\n",
      " [ 3.14692295e-01]\n",
      " [-6.37876816e-03]\n",
      " [ 3.57721943e-02]\n",
      " [ 1.41317555e-02]\n",
      " [-1.03194760e+00]\n",
      " [-3.63928712e-02]\n",
      " [ 5.07369262e-01]\n",
      " [-3.61951549e-01]\n",
      " [ 4.51190838e-01]\n",
      " [-1.08857178e-01]\n",
      " [-3.09652278e-01]\n",
      " [ 2.21163468e-01]\n",
      " [-7.55366414e-02]\n",
      " [-3.71034524e-01]\n",
      " [ 6.44307438e-01]\n",
      " [-1.16182383e-01]\n",
      " [ 1.55188141e-01]\n",
      " [-2.65133315e-01]\n",
      " [-1.80831914e-02]\n",
      " [ 2.78164571e-03]\n",
      " [-1.09021855e+00]\n",
      " [ 1.40483944e-01]\n",
      " [-1.98036461e-01]\n",
      " [ 2.47131415e-01]\n",
      " [-9.85062970e-01]\n",
      " [ 1.82448349e-01]\n",
      " [-4.45936399e-02]\n",
      " [-8.80224442e-01]\n",
      " [ 8.51014545e-02]\n",
      " [ 3.79698872e-02]\n",
      " [-4.99201125e-01]\n",
      " [-3.39445751e-02]\n",
      " [ 6.56062991e-01]\n",
      " [ 1.17035864e-01]\n",
      " [-4.60724149e-01]\n",
      " [ 7.31087065e-01]\n",
      " [-1.39063430e-01]\n",
      " [ 1.81208202e-01]\n",
      " [ 5.99665040e-01]\n",
      " [-4.70313585e-01]\n",
      " [ 3.46326923e-01]\n",
      " [ 4.89222268e-01]\n",
      " [ 3.84624021e-01]\n",
      " [ 3.39615269e-01]\n",
      " [-4.64103681e-01]\n",
      " [-2.86064207e-01]\n",
      " [-7.06135813e-02]\n",
      " [-1.52059731e-01]\n",
      " [ 3.39900697e-01]\n",
      " [ 1.02921046e-01]\n",
      " [-7.79282208e-01]\n",
      " [-8.41421597e-02]\n",
      " [ 1.41750971e-01]\n",
      " [-4.59267405e-01]\n",
      " [ 7.72793481e-01]\n",
      " [-8.95771733e-02]\n",
      " [-1.10278276e-02]\n",
      " [-3.46625585e-01]\n",
      " [ 4.86867406e-01]\n",
      " [ 5.01930166e-01]\n",
      " [ 3.47539937e-01]\n",
      " [-1.28217959e-01]\n",
      " [-6.46387250e-01]\n",
      " [ 6.39346579e-02]\n",
      " [-1.54177620e-01]\n",
      " [-1.30772405e-01]\n",
      " [ 1.61222532e-01]\n",
      " [ 1.42714368e-01]\n",
      " [ 2.61852687e-02]\n",
      " [ 2.99496906e-01]\n",
      " [-1.14192296e-01]\n",
      " [-7.44054381e-02]\n",
      " [ 3.61371990e-01]\n",
      " [-6.22519644e-01]\n",
      " [-1.42941257e-01]\n",
      " [ 2.00751781e-01]\n",
      " [-9.47084162e-02]\n",
      " [-5.95567445e-01]\n",
      " [ 3.89832513e-01]\n",
      " [-2.50390338e-01]\n",
      " [-1.39236134e-01]\n",
      " [-2.01327678e-01]\n",
      " [ 3.95343672e-01]\n",
      " [-5.67261498e-01]\n",
      " [ 1.50566153e-01]\n",
      " [-2.26056906e-01]\n",
      " [ 8.70591703e-02]\n",
      " [ 5.18856534e-02]\n",
      " [-2.86245497e-01]\n",
      " [ 1.35036017e-03]\n",
      " [-4.69838729e-01]\n",
      " [ 2.79179208e-01]\n",
      " [ 1.30075881e-01]\n",
      " [ 2.69777683e-01]\n",
      " [-5.63287839e-01]\n",
      " [ 5.62302630e-02]\n",
      " [-3.28763334e-01]\n",
      " [-2.23961213e-01]\n",
      " [-1.50770898e-01]\n",
      " [ 2.00771507e-01]\n",
      " [-2.94484277e-01]\n",
      " [-2.89186909e-01]\n",
      " [ 6.56945912e-01]\n",
      " [ 2.35401888e-01]\n",
      " [-2.14921931e-01]\n",
      " [ 5.13916096e-01]\n",
      " [ 2.68767061e-01]\n",
      " [ 2.33619308e-01]\n",
      " [-3.87701918e-01]\n",
      " [ 3.77047651e-01]\n",
      " [ 1.80517915e-01]\n",
      " [ 1.80105049e-01]\n",
      " [-4.47851876e-01]\n",
      " [ 9.47125812e-03]\n",
      " [ 2.47767101e-01]\n",
      " [-8.49820535e-01]\n",
      " [-1.45676958e-01]\n",
      " [-4.21794361e-01]\n",
      " [-4.70756816e-01]\n",
      " [ 4.10248873e-01]\n",
      " [-1.21209607e-02]\n",
      " [-4.48802318e-01]\n",
      " [-2.24866488e-01]\n",
      " [ 6.75531025e-04]\n",
      " [-5.74865333e-02]\n",
      " [ 1.40537911e-01]\n",
      " [-2.45560695e-02]\n",
      " [-3.52220876e-01]\n",
      " [ 8.75140147e-02]\n",
      " [ 4.82237153e-01]\n",
      " [ 1.24931414e-01]\n",
      " [-4.57352539e-01]\n",
      " [-2.39726655e-01]\n",
      " [ 3.56371302e-01]\n",
      " [-7.86655587e-03]\n",
      " [-2.95541874e-01]\n",
      " [-1.65370114e-01]\n",
      " [ 4.44903484e-01]\n",
      " [ 1.00318507e-01]\n",
      " [ 1.61567564e-02]\n",
      " [ 4.81695694e-02]\n",
      " [-5.19663129e-03]\n",
      " [-2.10197580e-01]\n",
      " [ 1.85474320e-01]\n",
      " [ 1.62119973e-01]\n",
      " [ 3.80853315e-01]\n",
      " [ 3.57998806e-01]\n",
      " [-5.21555402e-01]\n",
      " [-1.89306812e-01]\n",
      " [-4.13168656e-01]\n",
      " [-2.36699619e-01]\n",
      " [-1.68437135e-01]\n",
      " [ 5.43230648e-02]\n",
      " [ 3.48652254e-01]\n",
      " [-1.08794921e-01]\n",
      " [-3.40354489e-01]\n",
      " [-9.80998519e-02]\n",
      " [ 4.28125151e-01]\n",
      " [ 4.83603246e-01]\n",
      " [-2.92066753e-01]\n",
      " [ 8.76443817e-01]\n",
      " [-1.05546566e-01]\n",
      " [-3.60470946e-01]\n",
      " [-1.62395639e-01]\n",
      " [-3.50486405e-01]\n",
      " [-4.29574249e-02]\n",
      " [-5.16683828e-03]\n",
      " [ 7.04744561e-02]\n",
      " [-2.71839838e-01]\n",
      " [-2.83670044e-02]\n",
      " [-5.72040072e-02]\n",
      " [ 1.86420438e-01]\n",
      " [ 2.13876198e-01]\n",
      " [-4.80473688e-01]\n",
      " [ 3.84366401e-01]\n",
      " [ 3.60458398e-02]\n",
      " [-4.10821320e-01]\n",
      " [-3.19775201e-01]\n",
      " [ 5.41220258e-01]\n",
      " [-4.11812120e-01]\n",
      " [-7.07373665e-01]\n",
      " [ 7.07080372e-02]\n",
      " [ 1.36803176e-01]\n",
      " [-5.32305871e-01]\n",
      " [-1.58361340e-01]\n",
      " [-5.92652745e-01]\n",
      " [-5.33868835e-01]\n",
      " [-6.62447867e-02]\n",
      " [ 2.48856194e-01]\n",
      " [-8.91900077e-02]\n",
      " [ 8.62777201e-01]\n",
      " [ 2.66100530e-01]\n",
      " [-7.14309328e-02]\n",
      " [ 5.68036604e-01]\n",
      " [ 4.59112371e-01]\n",
      " [ 4.56307308e-01]\n",
      " [-4.27017021e-02]\n",
      " [ 4.64372352e-01]\n",
      " [-5.36513206e-03]\n",
      " [ 1.53091202e-01]\n",
      " [-4.95573736e-01]\n",
      " [-8.26680570e-02]\n",
      " [ 6.06380248e-01]\n",
      " [-6.70576822e-02]\n",
      " [-1.36046316e-01]\n",
      " [ 6.15672414e-02]\n",
      " [ 1.00932700e-01]\n",
      " [ 1.57086779e-01]\n",
      " [ 6.13154568e-01]\n",
      " [ 1.17778432e-01]\n",
      " [-1.31840498e-01]\n",
      " [-8.65177056e-01]\n",
      " [-2.51235081e-01]\n",
      " [ 1.58613914e-01]\n",
      " [ 3.00075679e-01]\n",
      " [ 1.38222003e-01]\n",
      " [ 3.22213789e-01]\n",
      " [-1.45384487e-01]\n",
      " [ 1.52445303e-01]\n",
      " [-2.01601844e-02]\n",
      " [ 1.48028106e-01]\n",
      " [ 3.90998946e-01]\n",
      " [-3.79953422e-01]\n",
      " [-3.48396076e-01]\n",
      " [ 9.53218658e-02]\n",
      " [ 3.25302709e-01]\n",
      " [-1.18009161e-01]\n",
      " [ 9.30146892e-03]\n",
      " [ 3.16684308e-01]\n",
      " [ 3.26472978e-01]\n",
      " [-2.79603943e-01]\n",
      " [ 2.16905057e-01]\n",
      " [-3.90755633e-02]\n",
      " [-5.20087254e-01]\n",
      " [-2.09884855e-01]\n",
      " [ 1.65923226e-01]\n",
      " [ 4.82235178e-01]\n",
      " [-9.72503050e-02]\n",
      " [ 2.79124110e-01]\n",
      " [-6.89751968e-02]\n",
      " [-6.79729170e-01]\n",
      " [ 2.32558641e-02]\n",
      " [-2.20847299e-01]\n",
      " [ 7.60908725e-02]\n",
      " [ 2.64826443e-01]\n",
      " [-1.62022949e-01]\n",
      " [ 1.09036569e-01]\n",
      " [ 3.03343125e-01]\n",
      " [-1.60051490e-01]\n",
      " [-2.53758404e-01]\n",
      " [ 4.81055731e-01]\n",
      " [ 1.36521453e-02]\n",
      " [-8.73299979e-02]\n",
      " [-1.10479600e-01]\n",
      " [ 1.47965271e-02]\n",
      " [-9.35948240e-01]\n",
      " [ 1.63114066e-01]\n",
      " [ 7.06095642e-01]\n",
      " [ 5.15013351e-02]\n",
      " [ 9.14893797e-04]\n",
      " [-3.14586406e-01]\n",
      " [-4.04070311e-01]\n",
      " [-4.03626178e-01]\n",
      " [ 6.92500298e-01]\n",
      " [-3.35364812e-01]\n",
      " [-3.39969803e-01]\n",
      " [-1.14262714e-01]\n",
      " [ 1.21278223e-01]\n",
      " [-4.08418448e-01]\n",
      " [ 1.58331527e-01]\n",
      " [ 1.11302771e-01]\n",
      " [ 4.77263071e-02]\n",
      " [-6.00150725e-01]\n",
      " [ 6.51415544e-02]\n",
      " [ 2.33280592e-01]\n",
      " [ 6.04800270e-01]\n",
      " [ 3.32565366e-01]\n",
      " [ 7.72005448e-02]\n",
      " [ 2.43383801e-02]\n",
      " [-9.88181853e-02]\n",
      " [ 3.78371150e-01]\n",
      " [-4.22041625e-01]\n",
      " [-1.46374955e-01]\n",
      " [-8.59745358e-02]\n",
      " [-1.19605023e-01]\n",
      " [-2.08852139e-02]\n",
      " [ 5.85533206e-02]\n",
      " [ 9.05670788e-02]\n",
      " [-1.96757628e-01]\n",
      " [ 3.56026066e-02]\n",
      " [-8.80625503e-02]\n",
      " [-9.93514901e-02]\n",
      " [-3.28487513e-01]\n",
      " [ 5.28228603e-01]\n",
      " [-3.07161835e-01]\n",
      " [-1.70332656e-02]\n",
      " [ 2.09444206e-01]\n",
      " [-7.37295626e-01]\n",
      " [-1.79936902e-01]\n",
      " [-2.08843002e-01]\n",
      " [-3.13803501e-01]\n",
      " [-5.43088958e-01]\n",
      " [-1.82396661e-01]\n",
      " [ 3.87972984e-01]\n",
      " [ 2.12647066e-01]\n",
      " [-1.24351356e-01]\n",
      " [ 9.20448431e-02]\n",
      " [ 1.62717135e-01]\n",
      " [ 6.33578072e-02]\n",
      " [ 3.47349130e-01]\n",
      " [-1.87101242e-01]\n",
      " [-5.35347614e-01]\n",
      " [-2.14613597e-01]\n",
      " [-1.15950870e-01]\n",
      " [-4.71298016e-01]\n",
      " [ 3.26461049e-01]\n",
      " [ 4.78024599e-01]\n",
      " [ 2.03337904e-02]\n",
      " [-5.87435142e-01]\n",
      " [-4.94025067e-02]\n",
      " [-9.94971099e-02]\n",
      " [-1.66832122e-01]\n",
      " [ 1.48252810e-01]\n",
      " [-3.36913910e-01]\n",
      " [-2.18728818e-01]\n",
      " [-3.17350529e-01]\n",
      " [ 2.42731831e-01]\n",
      " [ 3.32988190e-01]\n",
      " [-1.21163255e-01]\n",
      " [ 3.44173017e-01]\n",
      " [ 2.60635806e-01]\n",
      " [ 9.14279342e-02]\n",
      " [ 2.81869715e-01]\n",
      " [ 1.25171982e-01]\n",
      " [-5.65420001e-01]\n",
      " [ 1.60601048e-01]\n",
      " [ 2.28877220e-01]\n",
      " [-1.55642284e-01]\n",
      " [-2.94882916e-01]\n",
      " [ 8.30418356e-01]\n",
      " [ 4.86977705e-01]\n",
      " [-5.28225190e-01]\n",
      " [-5.55809328e-01]\n",
      " [-1.06756126e-01]\n",
      " [ 1.67183344e-01]\n",
      " [-1.61863169e-01]\n",
      " [ 4.03107456e-01]\n",
      " [-4.10361837e-01]\n",
      " [ 5.43039442e-01]\n",
      " [ 9.55214149e-02]\n",
      " [-2.62904291e-01]\n",
      " [ 2.44024121e-01]\n",
      " [-4.98088088e-02]\n",
      " [-3.38208278e-01]\n",
      " [ 1.03427125e-01]\n",
      " [ 6.80446677e-01]\n",
      " [-2.14374860e-01]\n",
      " [-8.14683617e-01]\n",
      " [-3.67053313e-01]\n",
      " [-1.80344720e-01]\n",
      " [ 2.79891301e-01]\n",
      " [ 4.10171924e-01]\n",
      " [ 4.34687763e-01]\n",
      " [ 1.36677498e-01]\n",
      " [-1.04074571e-01]\n",
      " [-4.55688936e-02]\n",
      " [-3.82539679e-01]\n",
      " [-4.25440589e-01]\n",
      " [ 8.77660813e-01]\n",
      " [ 5.58289096e-01]\n",
      " [-8.00714099e-01]\n",
      " [-3.20428389e-01]\n",
      " [-1.46059748e-01]\n",
      " [ 9.03035920e-01]\n",
      " [ 5.63613677e-01]\n",
      " [ 2.96063700e-01]\n",
      " [-3.39030770e-01]\n",
      " [ 4.37362208e-01]\n",
      " [-1.94243642e-01]\n",
      " [-1.11315348e-01]\n",
      " [ 2.66435578e-01]\n",
      " [ 7.45367077e-01]\n",
      " [-2.17344614e-01]\n",
      " [ 4.89332526e-01]\n",
      " [ 1.02000490e-02]\n",
      " [-1.16473930e-01]\n",
      " [-1.12839934e-01]\n",
      " [ 2.67856067e-01]\n",
      " [-3.45316700e-01]\n",
      " [ 7.51541538e-02]\n",
      " [ 5.26245276e-02]\n",
      " [ 4.17723021e-02]\n",
      " [ 2.03705984e-02]\n",
      " [-6.11748309e-01]\n",
      " [ 1.78637792e-01]\n",
      " [ 1.10750404e-01]\n",
      " [ 2.16028094e-01]\n",
      " [ 4.59871551e-01]\n",
      " [ 2.73225962e-01]\n",
      " [-7.75557083e-01]\n",
      " [-4.68738675e-01]\n",
      " [-4.52036516e-01]\n",
      " [-2.05581754e-01]\n",
      " [ 1.12298971e-02]\n",
      " [-1.17627258e-01]\n",
      " [ 1.88099471e-01]\n",
      " [-2.99370298e-01]\n",
      " [-4.90555511e-01]\n",
      " [-2.38549338e-01]\n",
      " [-5.74938661e-02]\n",
      " [ 3.78117410e-01]\n",
      " [-2.95390620e-01]\n",
      " [ 8.26084128e-01]\n",
      " [-2.54383665e-01]\n",
      " [-1.29846682e-01]\n",
      " [ 1.60764988e-01]\n",
      " [ 4.53011003e-01]\n",
      " [-3.05149437e-01]\n",
      " [-1.90142441e-01]\n",
      " [ 4.03086603e-01]\n",
      " [ 3.02061654e-01]\n",
      " [ 1.39162969e-01]\n",
      " [-4.18242039e-01]\n",
      " [-6.49157173e-01]\n",
      " [ 2.37645770e-01]\n",
      " [-2.73019542e-01]\n",
      " [-2.49164971e-01]\n",
      " [-2.07620234e-01]\n",
      " [-3.01992327e-01]\n",
      " [-1.96875159e-01]\n",
      " [-2.62274556e-01]\n",
      " [ 2.38194079e-01]\n",
      " [-2.89010997e-01]\n",
      " [-1.65345522e-02]\n",
      " [-2.40252435e-01]\n",
      " [ 2.47391016e-01]\n",
      " [ 2.12914767e-01]\n",
      " [-3.14597488e-01]\n",
      " [ 1.98231992e-01]\n",
      " [ 1.31882314e-01]\n",
      " [-8.29125883e-02]\n",
      " [ 6.15423350e-01]]\n"
     ]
    }
   ],
   "source": [
    "print(theta)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "a2f8a235",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x21b6bdb3588>]"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhpklEQVR4nO3de3Cd9X3n8ffn6OL7XbIRvmBjzMXcbHBsJ00yadJwy86azKYJtNswlMTdGeimO21nod0dUlo6m9006aZD6JLBDdAklDbN4KZ0iUNJINuAEWDAF2wLG2ML25ItW5YvsiWd7/5xHsknjmTd9Ujn+bxmzpzn/J7nOef707G/evT7/Z7fTxGBmZllQy7tAMzMbOQ46ZuZZYiTvplZhjjpm5lliJO+mVmGlKcdwPlUVVXFwoUL0w7DzGxMefXVVw9FRHV3+0Z10l+4cCG1tbVph2FmNqZI2tPTPjfvmJlliJO+mVmG9Jr0JY2XtFHSG5K2SPqTpPzbknZL2pQ8liXlkvQNSXWS3pR0XdF73SFpZ/K4Y9hqZWZm3epLm/5p4OMRcVxSBfAzSf+S7PvDiPiHc46/GViSPFYBDwOrJM0E7gdWAAG8Kml9RBwZioqYmVnver3Sj4LjycuK5HG+CXvWAI8n570ETJdUA9wIbIiIpiTRbwBuGlz4ZmbWH31q05dUJmkT0EAhcb+c7HowacL5uqRxSdlcYG/R6fuSsp7Kz/2stZJqJdU2Njb2rzZmZnZefUr6EdEREcuAecBKSVcB9wGXAx8AZgL/dSgCiohHImJFRKyoru52mKmZmQ1Qv8bpR8RRSc8DN0XEV5Pi05L+BviD5HU9ML/otHlJWT3wsXPKfzKAmM26FRF05IP2fNDWkaetI2jvyNOWT547gvZ8nrb2wnM+CufkAzry0bWdjzj7yHe+Ljo2Oo8t3n/23EIsRW2gnWVFL6O7sqJ6/GK9IDj7vr983i/uK36P8x0/UIOdjX1Qpw/ywy+unsyty3+pgSFTek36kqqBtiThTwA+CXxFUk1E7Jck4FZgc3LKeuAeSU9S6MhtTo57FvhzSTOS426g8NeCZURE0NqWp6W1jWOt7bS0ttHS2p48ku3T7Zw6005rW55TbR20tnXQ2pZPnjtobe/g1JlC2en2Ds60538hyVvpkwZ2XgTkBGuWXYgG+iYloC9X+jXAY5LKKDQHPRURP5T0r8kvBAGbgP+UHP8McAtQB5wE7gSIiCZJfwq8khz3QEQ0DVlNLDVtHXkONLfy/tFTvN98isaW0xw+foZDx8/QdOI0h0+cSV6f5nR7/rzvJcH48jImVJYxvjzH+Iqy5FHYnjqhggkVZYxLXleW5SjPiYryHBU5UV6Wo7xMVORyVJQVXleUifJcUl6WoyI5J5cTOUFOQoIynS2TRE6iLNmXk8jlkufO18lDglzu7LFd6USg5FVnjlFXPVW03bnv7Ml9Pp5fTILF+7p7j87Xg815YzFp/uWPd/CXP95JxODrP5b1mvQj4k1geTflH+/h+ADu7mHfOmBdP2O0UaC1rYNdjSfYdeg47zQUnvcdOcX7R09x8Fgr+XMusivLc1RNqmTW5HHMmlzJktlTqJpcyfSJlUwZX86U8eVMHV+RbFd0lU2qLCeXy/D/SBs2uSTTZ/3vwVE9946lo+nEGTbXN/NWfTNb3m9mc/0x9h452dWcKsGF0yYwf+YEPrS4irnTx3Ph9AnMnTGBmmkTmDN1HJPHlY/Jq0ErXZ3/GvMRlJHdf5tO+sb+5lP8/J3DvLTrMC/tauK9ppNd+xbMnMhVc6fyH66bx8XVk1hcPZlFVZOYUFmWYsRm/df5F2TWlwV30s+gM+15fr7rMBu2HuDFnYfYc7iQ5KdNqGDVopn85qoFXD13GldeOI1pEytSjtZsaOUznvWd9DOirSPPT7Y38vSmen66vZGW0+1MrCzjQ4ur+PwHF/LBi2dx+QVT3J5uJaurTT/bOd9Jv9Rt23+Mv6/dx9Ob6jl84gyzJlXyqWtquOHKOXxocRXjK9xMY9nQeT0z2PsUxjon/RKUzwc/2dHAt17Yzc93HaaiTPzaFXP4zPXz+Oil1VSUeUZty57OcQXnjjTLGif9EtKRD37wej0P/6SOdxpPcMHU8dx38+V8dsV8ZkyqTDs8s1Sdbd7JdtZ30i8BEcHz2xv4yr9sZ/vBFpbWTOUvP7eMT11T46t6s3P4St/GtK3vH+NP/mkLL+9uYuGsiTz0G9dxy9UXeIy82Tl8pV/gpD9GnTrTwf9+biffenEX0ydU8MCaK7l95QJf2Zv1oKsjN9s530l/LHpj71G+9OTrvHv4JJ9bMZ/7brmc6RPdZm92Pp1//Xqcvo0Z+Xzw1y+8w9d+tIPZU8bx3S+u4kOLq9IOy2xMODtkM9uc9MeIY61t/JcnN/Hc2w186uoa/vzTV/tuWbP+8JU+4KQ/JuxqPM4XH69lz+GTPLDmSn5r9UXuqDXrp66bzbOd8530R7tNe49y599sRBJP3LWKDy6elXZIZmNS5/oDHrJpo9bPdh5i7RO1VE0exxN3reSiWZPSDslszMp13ZGb7azvpD9KPf92A2ufqGVx9WQe/+2VzJ46Pu2QzMY0L6JS4KQ/Cv1b3SF+529f5bILpvCdu1a7w9ZsKHRe6We8fcd38owyr+45whcer2XhrIk8/turnPDNhkjOgx8AJ/1RZc/hE3zhsVeYPWUcf/uFVcz0JGlmQ6Z4ucQs6zXpSxovaaOkNyRtkfQnSfkiSS9LqpP0d5Iqk/Jxyeu6ZP/Cove6LynfLunGYavVGNR8qo3f/vYrBPA3d65k9hS34ZsNpVyS7TLeutOnK/3TwMcj4lpgGXCTpNXAV4CvR8QlwBHgruT4u4AjSfnXk+OQtBS4DbgSuAn4piSv4AG0d+S557uv8V7TSf76P17PoiqP0jEbap5wraDXpB8Fx5OXFckjgI8D/5CUPwbcmmyvSV6T7P+ECncSrQGejIjTEbEbqANWDkUlxrpvPLeTF3ce4s9uvYrVF3scvtlw8pV+H0gqk7QJaAA2AO8ARyOiPTlkHzA32Z4L7AVI9jcDs4rLuzmn+LPWSqqVVNvY2NjvCo01/6/uEH/1fB2fuX4en/vAgrTDMStZZztys531+5T0I6IjIpYB8yhcnV8+XAFFxCMRsSIiVlRXVw/Xx4wKDS2tfOnJTSyunswDa65MOxyzkublEgv6NXonIo4CzwMfBKZL6hznPw+oT7brgfkAyf5pwOHi8m7OyZyI4L7vv0VLaxsP/cZ1TKz0LRNmw+lsm37KgaSsL6N3qiVNT7YnAJ8EtlFI/p9JDrsDeDrZXp+8Jtn/r1HoOVkP3JaM7lkELAE2DlE9xpz1b7zPc2838Ic3XsZlF0xJOxyzkuchmwV9ubysAR5LRtrkgKci4oeStgJPSvoz4HXg0eT4R4EnJNUBTRRG7BARWyQ9BWwF2oG7I6JjaKszNjS2nOb+9VtYvmA6d/7KorTDMcsEL6JS0GvSj4g3geXdlO+im9E3EdEK/HoP7/Ug8GD/wywtX16/hZOnO/hfn7mGspzvEjQbCV4uscB35I6wF3c28s9v7ed3P34Jl8x2s47ZSJHb9AEn/RHV3pHngX/ayoKZE/niRy9OOxyzTDm7XGK2s76T/gj6zsvvsbPhOP/tU1cwvsI3I5uNJA/ZLHDSHyFHTpzhaxt28OFLqvjk0jlph2OWOe7ILXDSHyHf/EkdLa1t/Pd/t9Tr25qlwOP0C5z0R0DDsVYe//kePr18nsfkm6WkaxKGjGd9J/0R8M2fvEN7PvjPn7gk7VDMMsvLJRY46Q+z/c2n+O7G9/j16+d5YXOzFMnLJQJO+sPuoefriAju/lVf5ZulSV1DNrPNSX8YNbac5qlX9vGZ6+czf+bEtMMxyzTh0TvgpD+snvj5u7Tl83zxI55fxyxtnoahwEl/mJw608ETL+3hE5fP4eLqyWmHY5Z5uZyHbIKT/rD5/mv7OHKyzVf5ZqOEp1YucNIfBvl88OjPdnPNvGmsXDQz7XDMjKIJ11KOI21O+sPg+e0N7D50gi985GLffWs2Spydeyfbad9Jfxh8b+N7VE0ex81XXZB2KGaWODsNg5O+DaH9zaf417cb+OyKeVSU+cdrNlp49E6Bs9IQ+/vafeQDPveB+b0fbGYj5uw4/ZQDSZmT/hDqyAd/98pePnxJladcMBtluu7IzfilvpP+EHpxZyP1R09x20pf5ZuNNl5EpaDXpC9pvqTnJW2VtEXSl5LyL0uql7QpedxSdM59kuokbZd0Y1H5TUlZnaR7h6dK6fm7V/Yya1IlNyx1B67ZaOOO3ILyPhzTDvx+RLwmaQrwqqQNyb6vR8RXiw+WtBS4DbgSuBD4saRLk90PAZ8E9gGvSFofEVuHoiJpaz7VxnPbGvjN1QuoLPcfUGajjSdcK+g16UfEfmB/st0iaRsw9zynrAGejIjTwG5JdcDKZF9dROwCkPRkcmxJJP1ntxzgTEeeNcvO96Mxs7TkvFwi0M82fUkLgeXAy0nRPZLelLRO0oykbC6wt+i0fUlZT+XnfsZaSbWSahsbG/sTXqrWb3qfi2ZN5Np509IOxcy64SGbBX1O+pImA98Hfi8ijgEPA4uBZRT+EviLoQgoIh6JiBURsaK6unoo3nLYNbS08m/vHOLfX3uh78A1G7V8pQ99a9NHUgWFhP+diPhHgIg4WLT/W8APk5f1QPHwlXlJGecpH9P++c395APWLLsw7VDMrAc5X48BfRu9I+BRYFtEfK2ovKbosE8Dm5Pt9cBtksZJWgQsATYCrwBLJC2SVEmhs3f90FQjXU9vep8raqZyyWwvem42Wslt+kDfrvR/Bfgt4C1Jm5KyPwJul7SMQmf4u8DvAETEFklPUeigbQfujogOAEn3AM8CZcC6iNgyZDVJyd6mk2zae5R7b7487VDM7DxyXWvkphtH2voyeudnnJ2Kutgz5znnQeDBbsqfOd95Y9GzWw4AcMtVNb0caWZpynlqZcB35A7aj7Ye5PILprBgltfANRsLst6846Q/CE0nzlD7bhOfXDon7VDMrBe5nO/OAif9QXlu20HygaddMBsDvFxigZP+IGzYepCaaeO5au7UtEMxs16cvSM35UBS5qQ/QK1tHby48xC/dsUc35BlNgacnXsn21nfSX+AfrbzEKfaOrjhSrfnm40Fnlq5wEl/gH687SBTxpWzatGstEMxsz7IyZPvgJP+gEQEP93RyIeXVHkaZbMx4mxHbqphpM4ZawDeaTzO/uZWPnrp2JgQzsy8iEonJ/0B+OmOQwB8+JKqlCMxs75ym36Bk/4AvLizkYurJjF/pu/CNRsrPOFagZN+P51u7+ClXYfdtGM2xnhq5QIn/X6qffcIrW15PrLETTtmY4mv9Auc9PvphZ2NVJSJ1Rd7qKbZWOLlEguc9PvpxR2HuP6iGUwa16dFx8xslBCehgGc9Pvl0PHTbN1/jI8scXu+2VhzdvROtrO+k34/bNzdBMAHF7tpx2ysyXmOLMBJv1827m5iYmUZV8+dlnYoZtZPXVf6GW/fcdLvh5d2Heb6i2ZQUeYfm9lY4+USC3rNXpLmS3pe0lZJWyR9KSmfKWmDpJ3J84ykXJK+IalO0puSrit6rzuS43dKumP4qjX0jp48w/aDLaxcODPtUMxsALyISkFfLlnbgd+PiKXAauBuSUuBe4HnImIJ8FzyGuBmYEnyWAs8DIVfEsD9wCpgJXB/5y+KsWDj7iYiYJWHapqNSZ5ks6DXpB8R+yPitWS7BdgGzAXWAI8lhz0G3JpsrwEej4KXgOmSaoAbgQ0R0RQRR4ANwE1DWZnhtHF3E5XlOa6d7/Z8s7FInnAN6GebvqSFwHLgZWBOROxPdh0AOlcTmQvsLTptX1LWU/m5n7FWUq2k2sbGxv6EN6xe3t3E8vnTGVdelnYoZjZAOXmcfp+TvqTJwPeB34uIY8X7ovCrc0h+lBHxSESsiIgV1dWjYzz8sdY2trzf7KYdszEuJ3m5xL4cJKmCQsL/TkT8Y1J8MGm2IXluSMrrgflFp89LynoqH/VeffcI+YDVi9yJazaWyVf6fRq9I+BRYFtEfK1o13qgcwTOHcDTReWfT0bxrAaak2agZ4EbJM1IOnBvSMpGvY3vNlGeE8sXjJl+ZzPrhqTMd+T2ZQKZXwF+C3hL0qak7I+A/wE8JekuYA/w2WTfM8AtQB1wErgTICKaJP0p8Epy3AMR0TQUlRhur+05wpUXTmVCpdvzzcYy4Y7cXpN+RPyMs0Ncz/WJbo4P4O4e3msdsK4/AaatvSPPm/ua+dwH5vd+sJmNajnJ4/TTDmC0236whVNtHSxfMD3tUMxskHLyOH0n/V68/t5RAK5ze77ZmCfJHblpBzDavfbeEaomVzJvxoS0QzGzQZLwkM20AxjtNr13lGXzZ3TdzWdmY1ehIzftKNLlpH8eR06cYdehE1x30fS0QzGzIZDLKfOjd5z0z2PTvqMALJ/v9nyzUiB8c5aT/nm8vucIOcE18zzJmlkp8JBNJ/3zen3vUS67YKoXQTcrEZIy3o3rpN+jfD7YtPeox+eblRDJd+Q66fdgT9NJWlrbudZNO2YlwzdnOen36K36ZgCu8iLoZiVDuE3fSb8Hm+ubqSzLcemcKWmHYmZDxIuoOOn36K19zVxeM4WKMv+IzEqFp1Z20u9WRLD5/WY37ZiVGHfkOul3a8/hQifu1U76ZiUl5yGbTvrd6ezEddI3Ky2F5RKznfad9LvhTlyz0pRzm76Tfnfeqm/msgumUFnuH49ZKSnMvZPtrO+sdo6IYHO9O3HNSpF8c5aT/rneazrJsdZ2rpo7Ne1QzGyIFTpys531e036ktZJapC0uajsy5LqJW1KHrcU7btPUp2k7ZJuLCq/KSmrk3Tv0FdlaGyuPwa4E9esFEmQz6cdRbr6cqX/beCmbsq/HhHLksczAJKWArcBVybnfFNSmaQy4CHgZmApcHty7Kiz5f1mynPisgvciWtWanylD73OGRwRL0ha2Mf3WwM8GRGngd2S6oCVyb66iNgFIOnJ5Nit/Q95eL19oIXF1ZMZV16WdihmNgw8DcPA3SPpzaT5p3NpqbnA3qJj9iVlPZX/EklrJdVKqm1sbBxEeAOzbf8xLq/xVb5ZKSoM2cx21h9o0n8YWAwsA/YDfzFUAUXEIxGxIiJWVFdXD9Xb9snRk2fY39zKFTXuxDUrRbmcR+8MaEmoiDjYuS3pW8APk5f1wPyiQ+clZZynfNTYtr8FwEnfrER5auUBXulLqil6+Wmgc2TPeuA2SeMkLQKWABuBV4AlkhZJqqTQ2bt+4GEPj237CyN3rnDzjllJyomMd+P24Upf0veAjwFVkvYB9wMfk7SMws/vXeB3ACJii6SnKHTQtgN3R0RH8j73AM8CZcC6iNgy1JUZrG37jzFrUiXVk8elHYqZDQcp8x25fRm9c3s3xY+e5/gHgQe7KX8GeKZf0Y2wtw+0cEXNVCSlHYqZDYOcp1b2Hbmd2jvybD/Y4qYdsxIm3JHrpJ/YfegEZ9rz7sQ1K2E5uSPXST+xtasT10nfrFR5amUn/S5vH2ihokwsrp6cdihmNly8iIqTfqdt+4+xuHqy59A3K2Eesumk3+Xt/S1u2jErccLTMDjpA82n2jhwrNXLI5qVuFzOE6456QM7DxamX7jsArfnm5UyT7jmpA/AjoPHAVgy21f6ZqXOV/rGjoMtTKosY+70CWmHYmbDqLCISrY56VNI+pfMmUIu5+kXzEqZPA2Dkz4Umncune32fLNS5ztynfRpOnGGQ8dPe01cswzw3DtO+uxIRu4s8XBNs5InT63spN+Z9C9z0jcreZ5a2UmfHQdbmDK+nDlTvXCKWakrdOSmHUW6nPQPHufSOVO8cIpZBhSGbGY762c66UcEOw62ePoFs4yQfHNWppN+4/HTHD3ZxqVzPFzTLAvkIZu9J31J6yQ1SNpcVDZT0gZJO5PnGUm5JH1DUp2kNyVdV3TOHcnxOyXdMTzV6Z8dBwrTL/hK3ywbcvLcyn250v82cNM5ZfcCz0XEEuC55DXAzcCS5LEWeBgKvySA+4FVwErg/s5fFGnqHLnjpG+WDcKLqPSa9CPiBaDpnOI1wGPJ9mPArUXlj0fBS8B0STXAjcCGiGiKiCPABn75F8mI23GwhRkTK6iaXJl2KGY2AryIysDb9OdExP5k+wAwJ9meC+wtOm5fUtZTeao6O3E9cscsG9ymPwQduVG402HIfoqS1kqqlVTb2Ng4VG/7SyKCnclwTTPLBgny+bSjSNdAk/7BpNmG5LkhKa8H5hcdNy8p66n8l0TEIxGxIiJWVFdXDzC83u1vbqXldLtH7phliPBf9QNN+uuBzhE4dwBPF5V/PhnFsxpoTpqBngVukDQj6cC9ISlLjTtxzbInJ3fklvd2gKTvAR8DqiTtozAK538AT0m6C9gDfDY5/BngFqAOOAncCRARTZL+FHglOe6BiDi3c3hEOembZU9hucS0o0hXr0k/Im7vYdcnujk2gLt7eJ91wLp+RTeMdhw8TvWUccyY5JE7ZlkhX+ln947cwsgdt+ebZYm8XGI2k34+Xxi544XQzbLFyyVmNOnXHz3FqbYOr5ZlljE5T7iWzaS//UBnJ66bd8yypNCRm+2sn8mkv6PBSySaZVFh7p20o0hXNpP+gRZqpo1n6viKtEMxsxEkX+lnM+lvP3jcV/lmGeTlEjOY9Ns68rzTcJwr3Ilrljk5T7iWvaT/7qETnOnIe+SOWQYJT62cuaT/djJyx0nfLHtyOV/pZy7pbz/QQllOXDLbwzXNssZt+hlM+m8faGFR1STGlZelHYqZjTDhCdcymPSPcbmbdswyqbBcYrazfqaS/vHT7ew7cspJ3yyj5GkYspX0t3d14k5NORIzS4OHbGY06ftK3yyb5EVUspb0jzGpsoy50yekHYqZpaBzhdwsT8WQqaT/9oEWLrtgCrmcF0c2y6KcCv/3M5zzs5P0I4LtB1vcnm+WYUnOz3S7fmaSfkPLaY6ebHN7vlmG5bqSfrpxpGlQSV/Su5LekrRJUm1SNlPSBkk7k+cZSbkkfUNSnaQ3JV03FBXoqy3vNwNwRY2v9M2ySp3NOxkeqz8UV/q/GhHLImJF8vpe4LmIWAI8l7wGuBlYkjzWAg8PwWf32Vv7jiHBlRc66ZtlVWfzToZbd4aleWcN8Fiy/Rhwa1H541HwEjBdUs0wfH633qpv5uKqSUwaVz5SH2lmo4w7cgef9AP4kaRXJa1NyuZExP5k+wAwJ9meC+wtOndfUvYLJK2VVCuptrGxcZDhnbW5vpmr504bsvczs7Gnc9xeljtyB3vZ++GIqJc0G9gg6e3inRERkvr1042IR4BHAFasWDEk30xDSysHjrVylZO+WaZ1XemnHEeaBnWlHxH1yXMD8ANgJXCws9kmeW5IDq8H5hedPi8pG3ab6wuduNfMmz4SH2dmo5SHbA4i6UuaJGlK5zZwA7AZWA/ckRx2B/B0sr0e+Hwyimc10FzUDDSsNr13lJw7cc0yr2v0Tj7lQFI0mOadOcAPkh9iOfDdiPi/kl4BnpJ0F7AH+Gxy/DPALUAdcBK4cxCf3S8v727iygunuRPXLOM6x+lnecjmgLNgROwCru2m/DDwiW7KA7h7oJ83UKfbO9i09yi/ueqikf5oMxtlznbkphpGqkr+jtzN9c2cbs+zctGMtEMxs5R1zrvlCddK2MbdRwBYsXBmypGYWdp8pZ+BpP/CjkaWzJ5M1eRxaYdiZinr6sj1lX5pamhp5eXdh7n5qgvSDsXMRgGP0y/xpP/s5gPkAz51zYVph2Jmo4DH6Zdw0u/IB//wWj1LZk/mMk+nbGYUDdnMbs4v3aT/1R9t5429R/nCRxalHYqZjRJKunJ9pV9i6hqO839++g63r1zAZ1fM7/0EM8sET608+AnXRqVLZk/me19czfIFM7p6683M5KmVSzPpA6y6eFbaIZjZKJNzR25pNu+YmXXHQzad9M0sQzxk00nfzDLEbfpO+maWIZ3DOjwNg5lZBnS26XvCNTOzDPAiKk76ZpYhXR25GV4u0UnfzDKjqyPXV/pmZqXvbEduqmGkyknfzDIj5yGbI5/0Jd0kabukOkn3jvTnm1l2+easEU76ksqAh4CbgaXA7ZKWjmQMZpZdZ4dsZjfpj/SEayuBuojYBSDpSWANsHWE4zCzDOq80v/d773OhIqydIPpxeU1U/mr25cP+fuOdNKfC+wter0PWFV8gKS1wFqABQsWjFxkZlbyls2fzmeun8fJM+1ph9Kr+TMmDMv7jrqplSPiEeARgBUrVmT3bzAzG3LTJ1by1V+/Nu0wUjXSHbn1QPFSVvOSMjMzGwEjnfRfAZZIWiSpErgNWD/CMZiZZdaINu9ERLuke4BngTJgXURsGckYzMyybMTb9CPiGeCZkf5cMzPzHblmZpnipG9mliFO+mZmGeKkb2aWIRrNa0VKagT2DOItqoBDQxTOWJG1OmetvuA6Z8Vg6nxRRFR3t2NUJ/3BklQbESvSjmMkZa3OWasvuM5ZMVx1dvOOmVmGOOmbmWVIqSf9R9IOIAVZq3PW6guuc1YMS51Luk3fzMx+Ualf6ZuZWREnfTOzDCnJpJ+VxdclvSvpLUmbJNUmZTMlbZC0M3mekXacgyFpnaQGSZuLyrqtowq+kXzvb0q6Lr3IB66HOn9ZUn3yXW+SdEvRvvuSOm+XdGM6UQ+cpPmSnpe0VdIWSV9Kykv2ez5PnYf/e46IknpQmLL5HeBioBJ4A1iadlzDVNd3gapzyv4ncG+yfS/wlbTjHGQdPwpcB2zurY7ALcC/AAJWAy+nHf8Q1vnLwB90c+zS5N/4OGBR8m+/LO069LO+NcB1yfYUYEdSr5L9ns9T52H/nkvxSr9r8fWIOAN0Lr6eFWuAx5Ltx4Bb0wtl8CLiBaDpnOKe6rgGeDwKXgKmS6oZkUCHUA917ska4MmIOB0Ru4E6Cv8HxoyI2B8RryXbLcA2Cutpl+z3fJ4692TIvudSTPrdLb5+vh/mWBbAjyS9miwoDzAnIvYn2weAOemENqx6qmOpf/f3JM0Z64qa7UqqzpIWAsuBl8nI93xOnWGYv+dSTPpZ8uGIuA64Gbhb0keLd0bh78KSHpObhTomHgYWA8uA/cBfpBrNMJA0Gfg+8HsRcax4X6l+z93Uedi/51JM+plZfD0i6pPnBuAHFP7cO9j5p27y3JBehMOmpzqW7HcfEQcjoiMi8sC3OPunfUnUWVIFheT3nYj4x6S4pL/n7uo8Et9zKSb9TCy+LmmSpCmd28ANwGYKdb0jOewO4Ol0IhxWPdVxPfD5ZHTHaqC5qHlgTDunzfrTFL5rKNT5NknjJC0ClgAbRzq+wZAk4FFgW0R8rWhXyX7PPdV5RL7ntHuxh6ln/BYKveHvAH+cdjzDVMeLKfTmvwFs6awnMAt4DtgJ/BiYmXasg6zn9yj8mdtGoR3zrp7qSGE0x0PJ9/4WsCLt+Iewzk8kdXozSQA1Rcf/cVLn7cDNacc/gPp+mELTzZvApuRxSyl/z+ep87B/z56GwcwsQ0qxecfMzHrgpG9mliFO+mZmGeKkb2aWIU76ZmYZ4qRvZpYhTvpmZhny/wHMAFTvfniulwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "print(J_hist.size)\n",
    "x_axis = np.linspace(0,250,250)\n",
    "plt.plot(x_axis,J_hist[0:250])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97397ac5",
   "metadata": {},
   "source": [
    "# 2.\t对randn_data_regression_A的某一列乘上一个很大的数（比如100000），再调用梯度下降法计算线性回归模型的参数，你有什么发现？你觉得如何处理才有可能降低梯度下降法的迭代步数？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "id": "13943abf",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.         -1.66558438  0.61446305 ... -0.1933223  -0.02022444\n",
      "  -1.27698055]\n",
      " [ 1.          0.12533231  0.50774078 ...  0.07251873 -0.7443348\n",
      "   1.05665513]\n",
      " [ 1.          0.28767642  1.69242987 ...  0.94838331  0.16651915\n",
      "   0.3859412 ]\n",
      " ...\n",
      " [ 1.         -0.92190162  0.56896065 ... -1.47215039  0.31523604\n",
      "  -0.05946208]\n",
      " [ 1.         -2.17067449 -0.82171429 ... -1.40166474  2.57047553\n",
      "  -2.09209424]\n",
      " [ 1.         -0.05918782 -0.26560685 ...  1.03986473 -0.26376483\n",
      "  -1.36321714]]\n",
      "[[ 1.00000000e+00 -1.66558438e+00  6.14463049e+02 ... -1.93322299e-01\n",
      "  -2.02244400e-02 -1.27698055e+00]\n",
      " [ 1.00000000e+00  1.25332306e-01  5.07740785e+02 ...  7.25187290e-02\n",
      "  -7.44334802e-01  1.05665513e+00]\n",
      " [ 1.00000000e+00  2.87676420e-01  1.69242987e+03 ...  9.48383307e-01\n",
      "   1.66519145e-01  3.85941200e-01]\n",
      " ...\n",
      " [ 1.00000000e+00 -9.21901624e-01  5.68960646e+02 ... -1.47215039e+00\n",
      "   3.15236043e-01 -5.94620760e-02]\n",
      " [ 1.00000000e+00 -2.17067449e+00 -8.21714292e+02 ... -1.40166474e+00\n",
      "   2.57047553e+00 -2.09209424e+00]\n",
      " [ 1.00000000e+00 -5.91878250e-02 -2.65606851e+02 ...  1.03986473e+00\n",
      "  -2.63764834e-01 -1.36321714e+00]]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# 这里是求解 方法二\n",
    "\n",
    "# 第一部分: 设定传入参数\n",
    "n = rows\n",
    "m = cols\n",
    "X0 = theta\n",
    "\n",
    "A1 = np.array(A, copy=True)\n",
    "A1[:,2] = A1[:,2]*1000\n",
    "A1 = mat(A1)\n",
    "print(A)\n",
    "print(A1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "id": "53397524",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第0次迭代, cost = 243.768989\n",
      "第10次迭代, cost = 44.050490\n",
      "第20次迭代, cost = 10.552179\n",
      "第30次迭代, cost = 2.902350\n",
      "第40次迭代, cost = 0.855970\n",
      "第50次迭代, cost = 0.262709\n",
      "第60次迭代, cost = 0.082717\n",
      "第70次迭代, cost = 0.026514\n",
      "第80次迭代, cost = 0.008612\n",
      "第90次迭代, cost = 0.002826\n",
      "第100次迭代, cost = 0.000935\n",
      "迭代训练结束,迭代次数:100, 偏差值cost=0.000935\n"
     ]
    }
   ],
   "source": [
    "\n",
    "B = B\n",
    "theta = np.zeros((m,1))\n",
    "learning_rate = 0.01\n",
    "step = 2000\n",
    "tol = 0.001\n",
    "print_step = 10\n",
    "\n",
    "# 梯度下降方法\n",
    "(theta,J_hist) = gradientDescent(A,B,theta,learning_rate,step,tol,print_step)\n",
    "#print(theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33470c80",
   "metadata": {},
   "source": [
    "for line in locals()['In']:\n",
    "    print(line)\n",
    "    \n",
    "后面的都是草稿, 还有一定价值,没删除"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a02625d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#function [X,err,iter] = Gauss_S(A,X0,b,n,step,tol)\n",
    "            #迭代法求解好像必须是n*n型的正定?\n",
    "def jacobi(X,Y,tol,step):\n",
    "    n = X.shape[1] # 获取维度 python操作\n",
    "    m = X.shape[0]\n",
    "    print('m row =',m)\n",
    "    print('n col =',n)\n",
    "    theta0 = np.zeros(n)# 创建0向量\n",
    "    theta  = np.zeros(n)\n",
    "    print('step =',step)\n",
    "    for i in range(int(step)):\n",
    "        err = 0\n",
    "#         print('i=',i)\n",
    "        for j in range(m):\n",
    "            theta[j] = (Y[j]-X[j].dot(theta0)+X[j,j]*theta0[j])/X[j,j]\n",
    "        theta0 = theta\n",
    "        \n",
    "        err = np.linalg.norm(X.dot(theta)-Y)\n",
    "        if(err < tol):\n",
    "            print(\"迭代训练结束,迭代次数:\"+i+\", 偏差值:\"+err)\n",
    "    return (theta,err)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e40564e4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "(theta,err) = jacobi(A,B,tol,step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6f89bc2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f9ae286",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(A[1].size)\n",
    "print(A[1,1].size)\n",
    "print(A[1,:].size)\n",
    "print(X0[1].size)\n",
    "print(X0.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acfdb2a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "    print(A.dot(X0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42929e18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 草稿部分 \n",
    "#def jacobi(A,B,sigma,N):\n",
    "#     n = len(A)\n",
    "#     x0 = []\n",
    "#     x = []\n",
    "#     for i in range(0,n):\n",
    "#         x0.append(0)\n",
    "#         x.append(0)\n",
    "#     for k in range(1,N+1):\n",
    "#         R = 0\n",
    "#         for i in range(0,n):\n",
    "#             sum_ax = 0\n",
    "#             for j in range(0,n):\n",
    "#                 sum_ax = sum_ax + A[i][j] * x0[j]\n",
    "#             x[i] = x0[i] + ((B[i] - sum_ax)/A[i][i])\n",
    "#             if abs(x[i] - x0[i]) > R:\n",
    "#                 R = abs(x[i] - x0[i])\n",
    "#         if R <= sigma:\n",
    "#             print(\"精确度等于\",sigma,\"时，jacobi迭代法需要迭代\",k,\"次收敛\")\n",
    "#             return (x,k)\n",
    "#         for i in range(0,n):\n",
    "#             x0[i] = x[i]\n",
    "#     return (x,k)\n",
    "# def Jacobi(A,X0,B,n,step,tol):\n",
    "#     theta = X0\n",
    "#     for i in range(step):\n",
    "#         error=0\n",
    "#         for j in range(n):\n",
    "#             theta[j] = (B[j]-A[j,:]*X0+A[j,j]*X0[j])/A[j,j]\n",
    "# #             print(theta[j])\n",
    "# #             print(A[j,j])\n",
    "#         X0 = theta\n",
    "#         err = np.linalg.norm(A*theta-B-0)\n",
    "#         if(err < tol):\n",
    "#             print(\"迭代次数:\"+i+\", 偏差值:\"+err)\n",
    "# #AX = b Xθ=Y\n",
    "# # def = jacobi(A,B,tol,step):\n",
    "# #     n = A.ndim\n",
    "# #     theta0 = np.zeros(n)\n",
    "# #     theta  = np.zeros(n)\n",
    "# #     for i in range(step):\n",
    "# #         err = 0\n",
    "# #         print('i=',i)\n",
    "# #         for j in range(n):\n",
    "# #             theta[j] = (B[j]-A[j].dot(theta0)+A[j,j]*theta0[j])/A[j,j]\n",
    "\n",
    "# #         err = np.linalg.norm(A*theta-B-0)\n",
    "\n",
    "        \n",
    "            \n",
    "        \n",
    "# def jacobi(a,b,c=0.0001,d=30):\n",
    "#     x1=np.zeros(a.shape[1])\n",
    "#     x2=np.zeros(a.shape[1])\n",
    "#     k=0\n",
    "#     while k<d:\n",
    "#         k=k+1\n",
    "#         print('k=',k)\n",
    "#         for i in range(a.shape[1]):\n",
    "#             x2[i]=(-a[i].dot(x1)+b[i]+a[i,i]*x1[i])/a[i,i]\n",
    "#         if np.max(np.abs(x2-x1))<=c:\n",
    "#             print(\"x%d=\" % k,x2)\n",
    "#             print(np.max(np.abs(x2-x1)))\n",
    "#             break\n",
    "#         print(\"x%d=\" % k,x2)\n",
    "#         x1=x2.copy()\n",
    "#     return x2\n",
    "\n",
    "# def gauss_seidel(A,B,sigma,N):\n",
    "#     n = len(A)\n",
    "#     x0 = []\n",
    "#     x = []\n",
    "#     for i in range(0,n):\n",
    "#         x0.append(0)\n",
    "#         x.append(0)\n",
    "#     for k in range(1,N+1):\n",
    "#         R = 0\n",
    "#         for i in range(0,n):\n",
    "#             sum_ax = 0\n",
    "#             for j in range(0,n):\n",
    "#                 if j >= i:\n",
    "#                     sum_ax = sum_ax + A[i][j] * x0[j]\n",
    "#                 else:\n",
    "#                     sum_ax = sum_ax + A[i][j] * x[j]\n",
    "#             x[i] = x0[i] + ((B[i] - sum_ax)/A[i][i])\n",
    "#             if abs(x[i] - x0[i]) > R:\n",
    "#                 R = abs(x[i] - x0[i])\n",
    "#         if R <= sigma:\n",
    "#             print(\"精确度等于\",sigma,\"时，gauss_seidel迭代法需要迭代\",k,\"次收敛\")\n",
    "#             return (x,k)\n",
    "#         for i in range(0,n):\n",
    "#             x0[i] = x[i]\n",
    "#     return (x,k)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
